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ABSTRACT 

A  novel  transform  is  presented  which  maps  continuum 
functions  (such  as  probability  distributions)  into  disc¬ 
rete  sequences  and  permits  rapid  numerical  calculation 
of  convolutions,  multiple  convolutions,  and  Neumann 
expansions  for  Volterra  integral  equations.  The  trans¬ 
form  is  based  on  the  Laguerre  polynomials,  associated 
Laguerre  functions,  and  their  convolution  properties. 

Part  1  of  this  paper  deals  with  functions  having 
support  only  on  [0,5) .  The  resulting  unilateral 
Laguerre  transform  finds  applications  in  convolution  of 
such  functions,  inversion  of  Laplace  transform,  and  in 
solution  to  renewal  and  related  Volterra  integral 
equations.  —  '  ■*<  •  ‘  *  t  »  >  J 

Part  2  of  this  paper  deals  with  functions  having 
support  on  (-<»,<»)  via  a  bilateral  Laguerre  transform 
which  is  an  extension  of  the  unilateral  transform. 
Applications  of  this  technique  include  convolution  of 
such  functions  and  analysis  of  the  Lindley  process. 

Part  1  has  been  published  in  Applied  Mathematics 
and  Computation  and  part  2  has  been  submitted  for  pub¬ 
lication  in  that  journal. 
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LAGUERRE  TRANSFORMATION  AS  A  TOOL  FOR  THE  NUMERICAL  SOLUTION 
OF  INTEGRAL  EQUATIONS  OF  CONVOLUTION  TYPE 

J.  Keilson  and  W.  R.  Nunn 

ABSTRACT 

A  novel  transform  is  presented  which  maps  continuum  func¬ 
tions  (such  as  probability  densities)  into  discrete  sequences 
and  permits  rapid  numerical  calculation  of  convolutions,  multiple 
convolutions  and  Neumann  expansions  for  Volterra  integral  equa¬ 
tions.  The  transform  is  based  on  the  Laguerre  polynomials, 
associated  Laguerre  functions  and  their  simple  convolution  prop¬ 
erties.  A  second  transform  employs  Erlang  functions  as  elements 
of  the  basis.  The  limitations  and  advantages  of  the  two  trans¬ 
forms  are  discussed.  Numerical  inversion  of  Laplace  transforms 
relates  simply  to  the  Erlang  transform.  The  deconvolution  of 
two  functions,  i.e. ,  the  solution  of  a(t)  =  x(t)*b(t)  may  also 
be  obtained  quickly  in  this  way. 


INTRODUCTION 


One  often  encounters  in  applied  studies  integral  equations  [16] 
either  of  form 


f 

J0 


a(x  -  x')f(x,)dx'  =  b(x) 


(1) 


or  of  the  form 


f  ( 


x)  -  f  a(x  -  x')f(x')dx'  =  b(x) 


(2) 


where  a(x)  and  b(x)  are  specified  functions  and  f(x)  is  to  be 
found.  Equations  (1)  and  (2)  are  said  to  be  Vol terra  integral  eq¬ 
uations  of  convolution  type  of  the  first  and  second  kind  respective¬ 
ly.  The  Neumann  series  solution  of  (2)  has  the  form  [19] 


f(x)  =  b(x)  +  b(x) 


I  a(k)(x) 


(3) 


(k ) 

where  the  asterisk  denotes  convolution  and  a  (x)  is  the  k-fold 
convolution  of  a(x)  with  itself. 


The  entity  £  a(k) (x)  and  matrix  variants  associated  with 
0 

systems  of  integral  equations  of  convolution  type  arise  in  opera¬ 
tions  research  [6  ]r  engineering  [  7  ]f  and  biological  studies  [lo3* 
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Sometimes,  differential-integral  equations  give  rise  to  expressions 
such  as 


s  (t)  = 


CD 

I  e"XT 

k=0 


(U)k 

(k  +  1)  ! 


a  (k+1) 


(T) 


(4) 


which  describes  the  busy-period  density  for  certain  M|G|l  queueing 
systems  [is]. 

In  easy  cases  the  integral  equations  may  be  solved  analy¬ 
tically  via  Laplace  transformation,  and  full  answers  may  be  obtained 
when  the  Laplace  transforms  are  invertible.  More  often  than  not, 
such  transforms  cannot  be  inverted  and  expressions  such  as  (4)  are 
of  limited  value  when  they  cannot  be  evaluated  explictly.  The 
Laguerre  transformation  techniques  developed  in  this  paper  may 
then  be  of  value. 

The  deconvolution  problem  of  finding  f(x)  from  (1)  when 
a(x)  and  b(x)  are  known  numerically,  say,  is  particularly 
troublesome,  and  start-up  difficulties  described  below  may  make 
conventional  numerical  procedures  useless. 

The  Laguerre  transform  techniques  described  map  continuum 
functions  into  sequences,  and  map  the  continuum  convolution  opera¬ 
tion  into  lattice  convolution  of  these  sequences.  Such  discrete 
convolutions  are  well  matched  to  modern  computer  competence,  and 
the  inversion  mapping  back  to  the  continuum  is  direct. 


0-2 
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Laguerre  transformation  has  been  developed  as  a  tool  for  the 
solution  of  differential  equations  [ 12  3  -  The  applications  of 
interest  here  are  quite  different  and  new  tools  have  been  needed  to 
convert  the  underlying  simple  idea  into  a  flexible  working  proce¬ 
dure  adapted  to  computer  requirements. 

The  first  section  introduces  the  Laguerre  transform 

T  +  00 

T:  £  ( t)  -*•  (f  )„  in  a  form  convenient  for  our  needs.  One  has 
n  0 


0-3 
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where  4>(s)  is  the  Laplace  transform  of  f(r).  This  relationship 

permits  evaluation  of  for  many  important  f(r). 

Section  2  provides  simple  examples  of  the  transform,  and 

Section  3  discusses  the  structure  of  T^(u)  in  the  complex  u- 

plane.  Such  insight  into  structure  in  the  complex  plane  is  crucial 
to  many  of  our  algorithms  and  theorems. 

It  is  often  desirable  to  work  with  an  expansion  of  the  form 


f  (t) 


«  £ 
0 


.E+ 


n 


f"'  -T-  e 

n 


T 

n ! 


-T/2 


(8) 


whose  basis  functions  have  convolution  properties  of  comparable 
simplicity  to  the  Laguerre  functions.  The  mapping  T:  f(-r)  -*• 

will  be  called  an  Erlang  transform.  The  set  of  functions  so  repre¬ 
sentable  is  more  limited  than  that  for  the  Laguerre  basis.  They 
must  be  integral  functions  of  t  of  order  at  most  one,  so  that 
2 

f(t)=  e  T  for  example  is  excluded.  The  nature  of  the  Erlang 
transform  and  its  relationship  to  the  Laguerre  transform  is  des¬ 
cribed  in  Section  4. 

The  square  integrability  requirement  on  f(x)  for  the 
Laguerre  transform  is  no:  an  intrinsic  limitation.  An  exponential 
transformation  is  described  in  Section  5  which  avoids  such  diffi¬ 
culties. 


0-4 


Even  though  square  integrabiliry  suffices  in  principle  for 
the  Laguerre  transform,  practicability  of  the  method  requires 

DO 

fall  off  quickly  with  n  so  that 

0 

computer  time  is  not  excessive.  In  Section  6,  rapidity  of  disap¬ 
pearance  of  (f^+)  with  n  is  related  to  the  smoothness  and  con¬ 
centration  of  f(t).  In  particular,  it  is  shown  that  "rapidly 
decreasing"  f(-r)  for  which  (d/dt )  pf  ( x )  -*•  0  as  t  -*•  »  for  all 
non-negative  q,p  are  associated  with  rapidly  decreasing  (f  )  for 
K  Lt 

which  n  fn  -*■  0,  n  -+■  °°  for  all  positive  K. 

Algorithms  for  the  calculation  of  the  Laguerre  coefficients 
are  presented  in  Section  7.  Section  8  is  devoted  to  a  discussion 
of  the  deconvolution  problem. 

A  variety  of  numerical  examples  of  the  method  are  treated  in 
Section  9,  and  the  implementation  of  the  procedure  is  discussed. 

Section  10  describes  interpolation  methods  and  problems  when 
the  known  functions  are  known  only  numerically. 

A  final  section  deals  with  possible  generalizations  of  the 
method  to  special  families  of  functions. 


that  the  coefficients 


(0 


0-5 
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1.  The  Laguerre  transform. 

Let  f{x)  be  any  function  in  L2(o,<*>),  i.e.,  any  square 

integrable  measurable  function  on  (o,°°).  Then  f  (x)  may  be  ex¬ 
panded  in  terms  of  the  Laguerre  functions 


Vx)  =  e_35X  Ln(x)  ,  (1) 

where  Lr(x)  are  the  Laguerre  polynomials  having  the  Rodrigues 
formula 


1_ 

n! 


»x  fL 
dx 


n 


/  n  -  x  ■> 
lx  e  } 


(2) 


Classical  properties  of  the  Laguerre  polynomials  and  functions  are 
given  in  an  appendix.  The  Laguerre  functions  £n(x)  are  ortho¬ 
normal  on  (0 , “) ,  i.e.. 


f  l  (xU  (x)  dx 
/  m  n 


,0,  m^n-> 
ll,  m=n'  ' 


(3) 


and  provide  a  complete  basis  for  L2  in  the  associated  Hilbert 
space  metric.  For  f (x)  e  L2,  one  has  the  Fourier-Laguerre 
expansion 


f(x)  =  Z  f!*n(x>»  (4) 

_  n  n 


1-1 
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where,  from  (3) 


i-r 


f(x)  £  (x)  dx  . 
n 


(5) 


Equation  (5)  describes  a  mapping  of  the  function  f(x)  on 
the  continuum  (o,“)  into  a  sequence  (f^.1  ,  i.e.  into  a  function 


(€ ■ 


on  the  non-negative  integers.  This  sequence  will  be  called  the 
Laguerre  transform  of  f(-r).  It  should  be  noted  that  "Laguerre 
transform"  is  often  used  to  denote  the  mapping 

JF(x)e”X  L  (x)d>:  =  /  F(x)e~X//2  t  (x)  dx.  This  mapping  has 

o  n  Jo  n 

been  employed  as  a  tool  for  the  study  of  differential  equations  [  11] 
From  (5)  one  has  the  Parseval  relation 


?  (f*)2  =  L  f2(x)dx  *  (6) 

Two  useful  results  from  analysis  are  needed.  The  Laguerre  func¬ 
tions  have  the  generating  function 
00 

T\  l  (x)  un  =  (l-u)~'exp  {-—  x  (1+u) (1-u)  }  ;  o  <  u  <  1,  (7) 

n  /. 

o 


-8- 

and  Laplace  transform 

V'  *  J'  *'*\Md*  ‘ik  <f^>n  <8> 

From  (8)  one  may  verify  the  basic  convolution  property  of  the 
Laguerre  functions 

Z  (x)  *£  (x)  =  f*  l  (x-yU(y)dy  (9) 

n  m  JQ  n  m 

=  Vx)*£m+n  (x)  *  Wx)  "  Wl(x) 


This  convolution  property  underlies  the  Laguerre  transform  and  its 


utility.  For  any  f(x) 

in  L2  with 

Laguerre  transform  ^f* 

I 

we  define  a  pair  of  related  generating 

and  the  transform  sequence  (fw)  by 

no 

functions  T^(u)  and 

T*(u) 

(u)  =  t  . 

o 

0  <  J  U  j  <  1  ; 

(10a) 

tJ(u>  -  = 

o 

(1-U)  £  f^U11, 

o 

o  <  j  u  |  <  1  . 

(10b) 

The  reason  for  the  factor  (1-u)  in  (10)  will  soon  be  clear.  From 

+  i 

(6)  one  has  f  -*•  o,  as  n-*-®,  so  that  T^  (u)  is  regular  in  the 


1-3 
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interior  of  the  unit  circle  ju|<-  1  in  the  complex  u-plane.  From 
(10b)  one  has 


f*  -  fl  -  fl  i  , 


n  >  1 


f*  =  f  + 
o  o 


(11) 


so  that  the  sequence  (f  )  is  square  summable  as  well.  Employment 
of  (5)  in  (10)  gives 

.  00 

Tj(u)  =  (1-u)  £  ^ 

o 


/: 


f(x)£  (x)  dx  ,  o  <  u  <  1. 


n 


When  f(x)  is  integrable  on  (o,»)  ,  since  Un(x)  |^1  the  order 
of  summation  and  integration  may  be  interchanged.  Hence  from  (7) 

one  has 


T*(u)  =  f  f(x)exp  {-h  x  (1+u) (l-u)'1} 


dx 


a.e. , 


T*(u) 


=  $  (*5  (1  +  u)  (1-u)  *)  ,  o  <  u  <  1  , 


(12) 


/*> 

C  V 

e  f(x)dx  is  the  Laplace  transform  of  f (x) . 

u 

f*g  «  f  f  (x-y)  g  (y)  dy.  Since  (s)  =  $f(s)  4>  (s )  ,  one 


let 


has  from  (12) 


1-4 


r 


l 

I 


1 
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Tf*g<U>  “  TJ<U>  Tg<U>  »  U3) 

for  functions  f,g  that  are  both  integrable  and  square  integrable 
on  (o,®),  It  follows  from  (13)  moreover  that 


(f*g> 


# 

n 


T  f* 
o  n- 


m 


(14) 


The  value  of  the  Laguerre  transform  is  now  clear.  The  transformation 

£ 

via  (9)  and  (11)  maps  functions  f  (x) ,  g(x)  into  sequences  (f^) , 
(g*) ,  and  their  continuum  convolution  f(x)*g(x)  is  mapped  into 
a  lattice  convolution,  and  thence  back  onto  the  continuum  via 


4.  | 

f  *  £  f  and  the  representation  (4) .  The  advantage  of  the  procedure 
n  o  m 

for  numerical  convolution  of  two  functions  known  only  numerically 
is  moot.  But  its  value  for  iterative  convolutions  and  multiple 
convolutions  of  a  function  f (x)  with  itself,  and  weighted  sums  of 
such  multiple  convolutiorts  (i.e.  (0.4))  via  machine  computation  is 
apparent.  Some  of  the  applications  of  such  computation  in  statistics 
and  applied  probability  will  be  presented  subsequently.  The  trans¬ 
form  also  has  analytical  value,  as  we  will  also  attempt  to  show. 


l-5 
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2.  Some  Elementary  Laguerre  Transforms 
Exponential  and  Erlang  Functions 

The  simplest  and  most  basic  set  of  Laguerre  families  is  that 

9  T 

for  the  exponential  functions  f  { t)  =  e  .  Here,  since  <Hs)  = 
(Q+s)  \  we  have  from  (1.12) 


T 


# 

f 


(u) 


(1-u) 

(6  +  35)  -  (6-3s)u 


(1) 


and  hence 

The  Erlang  densities  for  scale  parameter  6  are  correspondingly 
simple,  since  they  arise  from  convolutions  of  the  exponential  den¬ 
sities.  They  are  given  in  the  table  of  transforms. 

Laguerre  Functions 

A  related  set  of  transforms  of  special  interest  is  that  for 
the  Laguerre  functions.'  Here  from  (1.10)  we  see  that 

L  (t)  t!(u)  »  (l-u)u11  (3) 

n  i 

in  particular 

*  <t>  =  e‘T/2  — *•  T*(u)  =  (1-u)  (4) 


2-1 
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so  that  for  its  N-fold  convolution 


-t/2 

Tn^TT  !  e 


tJ(u)  *  (l-u)N 


(5) 


For  this  family  of  functions  and  finite  mixtures  thereof,  the  trans- 
form  g.f.'s  T^(u)  are  polynomials  and  the  transform  sequences 
f*  terminate.  The  value  6=%  thus  provides  a  natural  scale  for 
computation  in  which  convergence  is  most  efficient.  This  will  be 
discussed  further  subsequently. 


2-2 


oM« 


3 


Structure  of  T^'(u)  in  the  complex  plane. 

When  f(t)  e  L2  (o ,«),  Equation  (1.6)  implies  that 

+  2  +  “  4-  n  J. 

(f^)  <  so  that  fn  -*■  0,  as  n  ®,  and  £  f  u  *  Tl(u)  is 

o 

regular  inside  the  unit  circle,  i.e.  for  0  <_  |u|  <  1.  Moreover 

4  # 

f  -*■  0  (and  is  square  summable)  so  that  (u)  is  regular  for 

0  £  | u ]  <  1,  as  may  also  be  seen  from  T*(u)  =  (l-u)T^(u).  It  may 

be  noted  that  ju|  =  1  maps  into  the  line  Re(s)  =  0,  i.e.  the 

imaginary  axis  in  the  s-plane,  and  that  the  point  u  =  1  corresponds 
to  s  = 

For  the  representation  f(x)  =  I  f+£  (t) ,  where  I Z  (t) |  <  1, 

n  n  n  — 

speed  of  convergence  is  promoted  by  a  geometric  decay  rate  for  f^, 

t  #  .  .  . 

and  hence  by  regularity  of  Tf (u) ,  Tf(u)  in  a  region  containing 

the  unit  circle  in  its  interior.  This  requires  (cf.  (1.12))  that 

<j>(s)  be  regular  at  s  *=  <*>  and  have  a  negative  abscissa  of 

convergence  [20]  •  Such  behavior  in  the  s-plane  of  4>(s)  also 

*  #  + 

assures  regularity  of  T^(u)  and  T^(u)  in  and  on  the  unit  circle. 
The  situation  is  summarized  by  the  following  theorem. 

Th"m  If  f(x)  eTT  e  L 1  (o,»)  for  some  y  >  o  and  if 

♦  (w*1)  =  f  f(x)e~T//w 
•'o 


dx  is  regular  at  w=o,  and  vanishes  there. 


-is- 


4.  The  Erlang  Family  of  Functions,  and  the  Erlang  Transform. 

The  class  of  functions  1^(0,°°)  is  too  broad  for  some  numerical 
work.  The  functions  one  deals  with  are  often  highly  smooth  and  such 
smoothness  permits  representation  by  sequences  (a^)  which  fall  off 
quickly  as  n  -*■  “.  An  alternative  to  the  Laguerre  building  blocks 
may  then  be  available  as  we  next  describe. 

T®  _  T  /  0 

Definition.  Let  e  (t)  =  -r  e  and  let  f  be  ^  family 

of  functions  (entire  functions  as  we  will  see) 


E  =  (a  (t)  :  a  ( t )  =  J2  a+e  ( t )  } 


where 


T^a*  <  °°,  some  6  >  0. 

m 


This  family  of  functions  will  be  called  the  Erlang  family. 

As  will  be  seen  condition  (lb)  permits  use  of  generating  func¬ 
tions  and  discrete  convolution  related  thereto.  We  note  that  E  is 
a  linear  space.  It  follows  from  (la),  and  (lb)  that 

CD 

a(t)eT//2  =  l  (a*/mi )  t111  is  a  power  series  with  infinite  radius  of 
0  ™ 

convergence  and  is  therefore  an  entire  function.  Hence  a(x)  is  also 
entire. 


4-1 
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Def inition.  Let  a(T)  e  E.  The  sequence  (a*)°°  will  be 
.  m  o 


Since  a(x)eT/^  =  £  a^Tm/m!  is  entire  and  hence  regular  at 

o 

t  =  0,  it  follows  from  the  Taylor  expansion  that 


Equation  (2)  provides  the  Erlang  expansion  of  any  function  a(x)  in 
the  Erlang  family  from  its  Taylor  expansion  about  x  =  0. 

The  Laplace  transform  of  a(x)  is  given  from  (1)  by 


a(s)  =  £a*  (s+}s)-m~1  (3) 

and  this  series  will  converge  absolutely  when  the  real  part  of  s  is 
sufficiently  large.  Let 


Ta+ (w) 

a 


„ .  m 
a*w  ; 
m 


T^*(w)  =  w 


5>*w 

m 

o 


(5a) 


The  transform  might  equally  be  called  a  Maclaurin  transform.  The 
name  Erlang  has  been  used  because  the  Laplace  transform  structure 
of  the  Gamma  densities  and  Erlang  densities  suggested  the  method. 
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i.e. 

Ta,<"’  -  a  (t£)  ' 

Hence  we  again  have,  for  a(-r),  b(r)  c  E 


Talb(w)  =  Ta#(w)  Tb#(w) 


(5b) 


(6) 


Again,  as  in  Section  1,  (6)  permits  study  and  evaluation  of  continuum 
convolution  via  lattice  convolution. 

E# 

The  validity  cf  the  transform  T  (w)  requires  that  this  be 

cl 

regular  in  w  at  w  =  0,  and  this  will  be  assured  when  a(s)  is 
regular  at  °°.  Indeed  such  regularity  provides  an  alternate  charac¬ 
terization  of  members  of  the  Erlang  family,  as  we  see  next. 

Theorem  4 . 1 

A  B 

,  .  r  <j=>  The  Laplace  transform  u(s)  of  a(t)  is 
a'  <’=>  regular  at  infinity  and  vanishes  there. 

Proof.  We  have  seen  that  A  «=>  B.  Suppose  conversely  that  B  is 

true.  Then  a(s)  regular  at  »  implies  a(s-h)  regular  at  <*>, 

i  .e. 

» 

ats-^)  =  f{  a(x)ei//2}  e  s‘dT  =  ^  C^d/s)™*1 


where  the  power  series  on  the  right  has  a  positive  radius  of 

oo 

convergence.  But  £  Cm  is  the  Laplace  transform  of 

o 
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w  . 

Y2  (C^/mDt  which  must  equal  a(x)eX  since  Laplace  transforms 
o 

are  one  to  one,  and  the  conclusion  follows. 

Closure  properties  of  the  space  E  are  also  of  interest,  and 
follow  at  once  from  Theorem  4.1. 

Theorem  4.2.  Let  a^(t)  ,  a2(t)  e  E.  Then 

(1)  c1a1(t)  +  c2a2(t)  e  E  ; 

(2)  a^tJe^1  c  E  ,  for  all  real  y  ; 

(3)  a^(.t)  *  a2(t)  €  E  (closure  under  convolution) 

We  observe  that  the  Laguerre  functions  Z (t)  and  finite 

linear  combinations  thereof  are  all  elements  of  E.  But  other 

elements  of  E  need  not  be  in  l2(o,®)  nor  even  go  to  zero  as 

2  -  3 

t  -*•  Thus  for  a(x)  *  x  ,  a(s)  =  2s  which  is  regular  at 

infinity.  Similarly  every  polynomial  in  x  is  in  the  space  E. 

It  may  be  seen  that  the  order  of  a(x)  is  at  most  equal  to 

1.  Hence  E  is  a  proper  subset  of  the  set  of  entire  functions  and 
_t2 

e  cannot  be  in  E. 

Note  that  the  space  E  is  not  closed  under  limits.  Thus 
aN(s)  *  (1  +  s/N)“N  converges  to  e‘s,  i.e.  a  sequence  of  Laplace 
transforms  regular  at  °°  converges  to  one  which  is  not. 
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1 


J- 
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A  subspace  of  E  of  special  interest  is  EL2  the  intersection 
of  E  and  L2(.o,^).  I£  a(\)  t  EL2,  it  will  have  both  an  Erlang 
transform  and  a  Laguerre  transform.  If  the  Laguerre  related  func¬ 
tion  T*(u)  of  (1.10b)  is  designated  by  (u)  =  a  ^  i we  have 

from  5b) 


Ta#(w)  =  T^#(l-w), 

a  a 


‘a 


(1-u) 


(7) 

(8) 


The  module  functions  (e  (t)}  and  {t  ( t ) }  of  the  two 

n  n 

transforms  are  easily  related.  From  (7), (8),  since  en^T)  = 

U  (t)] (n+1)  we  have  T^# (w)  »wn+1,  and  T^#(u)  =  (l-u)n+1  so 

n  n 

that 

n 

e  (r)  =  T  6  £  (x)  (9 

n  nr  r 

o 

where  6nr  =  {^)(-l)r.  Similarly 


VT) 


n 

-  z 

o 


:er(T; 


(10) 


That  (6  )  is  its  own  inverse  is  a  classical  result.  (See  for 

mn 

example,  Riordan  [15].  If  one  has  a(t)  c  EL2,  then 


4-5 


0.0 


(T>  *  ?aX(,)  *  Sai5em(T) 


(ID 


+  ,  « 


and  the  sequences  U  )  and  (a*)  are  interrelated  by 

no  mo  J 


a*  =  V  aT9 
m  J-'  n  nm 
n=m 


(12a) 


+  00 
i  =  a*  6 

n  J—1  m  mn 
m=n 


(12b) 


For  any  such  function  in  EL2  the  Erlang  transform  coefficients 
a*  obtained  from  (2)  provide  a  stepping  stone  to  the  Laguerre 

+  CO 

transform  (a  )  via  (12b). 
n  o 

The  Erlang  transform  method  for  evaluating  convolutions  is 
closely  related  to  the  method  of  phases  employed  by  Neuts  and  others 
[13,14]  for  the  numerical  evaluation  of  functions  arising  in  stochas¬ 
tic  models.  Indeed  the  method  of  phases  provides  an  intuitive  frame¬ 
work  for  understanding  the  Erlang  transform  method.  The  latter 
method  is  more  general  than  the  method  of  phases  in  several  respects. 

a)  It  is  applicable  to  convolutions  of  functions  of  mixed 
sign,  and  hence  to  the  response  functions  of  electrical  engineering 
arising  in  circuits  and  propagation,  for  example. 

b)  For  positive  functions  no  auxiliary  Markov  chains  models 
are  required.  Indeed,  certain  simple  probability  density  functions. 
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such  as  Ke-aT [1-cos  bij,  cannot  be  related  to  chain  models,  as  pointed 
out  by  D.  R.  Cox  [4].  Such  functions  present  no  difficulty  for  the 
Erlang  transform  method. 

c)  The  method  of  phases  requires,  in  principle,  integrable 
functions.  A  function  such  as  IQ(t)e  t,  where  Iq(t) 

Bessel  function,  is  unacceptable.  The  Erlang  method 
even  for  functions  which  grow  exponentially. 


is  the  modified 
may  be  employed 
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5.  Extension  of  the  LAguerre  Transform  to  Non-Inteqrable 

Functions  Via  Exponential  Transformation 

Our  applications  of  motivating  interest  center  about  the 

convolution  operation  d(t)  =  a(t)  *  b(t).  When  a(t)  and  b(t) 

are  not  in  I2 they  do  not  have  Laguerre  transforms.  We 

note  however  that  when  a(t)  and  b(t)  are  Laplace-transformable, 

one  has  6(s)  =  cx(s)$(s)  for  s  sufficiently  large.  Then 

a(s+8)B(s+8)  =  <$(s+8),  8  >  0,  i.e.  a0(t)  *  b0<t)  =  d0 (t) ,  where 

,  0 

(t)  *  a(t)e  ,  etc.  When  a0,  b0  e  they  may  be  convolved  via 
Laguerre  transform  methods  to  give  d0  and  one  may  then  find  d(t) 
from  d(t)  =  d0(t)e  .  Such  exponential  transformation  extends 
thereby  the  scope  of  Laguerre  transform  methods.  Note  that  the  same 
procedure  may  be  employed  for  the  integral  equations  of  convolution 
type  of  the  introduction.  For  example  the  equation 

h  (t)  =  b  (t)  +  a  ( t)  *  h  (t) 


becomes 


h0  (t)  =  b0(t)  +  a0  (t)  *  h0(t) 


under  such  transformation,  and  the  Laguerre  transform  method  may  be 
used  when  a0,bg  e  12  even  though  a,b  t  l7. 


5-: 
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6.  The  rate  of  convergence  of  the  Laguarre  coefficients 

It  is  well  known  that  a  Fourier  series  expansion  of  a  peri¬ 
odic  function  f(t)  =  eint  has  coefficients  a  which  fall 

n  n 

off  rapidly  with  n  when  f(t)  is  sufficiently  smooth  [53*  A 

•i.  $ 

similar  behavior  for  the  coefficients  f '  of  (1.4)  and  f  of 

n  n 

(1-11)  is  available  and  of  key  importance  to  the  methods  we  desire. 
We  first  observe  that  the  square  summability  of  (f^) ,  i.e., 

oo 

S'  (f+)^  <  »  implies  that  (f^)^-*-0  as  n  <*>.  Our  methods  re- 
n  n 

o 

"t 

quire  convolution  of  sequences  (f^) ,  (g^) .  As  remarked  in  Section 
3,  for  accuracy  of  calculation  we  would  like  the  coefficients  to 
fall  off  rapidly,  and  geometric  decay  of  the  coefficients  would  be 
desirable. 

It  may  be  noted  that  summability  £  | f ^ |  <  ®  implies  the 

o 

continuity  of  f(t)  in  (1.4).  The  argument  is  along  classical 
lines.  Since  |fn(t)  I  1  1/  one  has  |fn(t  +  h)  -  fn(t)  |  £  2. 

Hence  by  dominated  convergence 

00  + 

lim  f(t  +  h)  -  f(t)  =  lim  £  ^  Un(t  +  h)  -  *n(t)}  »  0. 

h-*-0  h+o  o 

and  f(t)  is  continuous. 

It  then  follows  that  when  f(t)  is  not  continuous  for  all 
t  in  that  (f m )  is  not  summ.able.  For  example  the 
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function  f (t)  =1#  0  £  t  £  1;  f (t)  *0,  1  <  t,  has  coefficients 

which  are  square  suromable  but  not  summable  and  hence  rather  slow  in 
convergence  to  zero.  Dominated  convergence  also  leads  to  the  con¬ 
clusion  that  f(x)  =  £f^£n(x)  goes  to  zero  as  t  -*•  »  when 

£  \ffn\  <  "• 

The  argument  can  be  generalized  to  the  following  theorem. 

Let  C+(R+)  be  the  linear  space  of  functions  f(t)  with  contin- 

(k) 

uous  k’th  derivatives  on  (o,“),  such  that  lim  f  (t)  *  0. 

x  -*•» 

“  + 

Theorem  2.1.  Let  f(t)  =  T  f_£_(T).  Then 
— . .  ■  "  *-*  mm 


V  mK|f+|  <  «  ->  f ( t )  £  Ck  (R,  )  »  k  =  1 , 2 , . . . . K.  The  proof 

m=0  m  1  + 

is  based  on  the  following  lemmas. 

Lemma  2.2  \l^]  (t)  i  <  \l^  (0)  ! 

m-1 

Proof.  The  standard  identity  L^(x)  =  "  £  L  (t)  implies 


that 


III  A  HI 

=  %£  (x)  +  T  lr{r)  =  V  6  If  t)  where  6  = 

m  m  —  r  r“g  mr  r  mm 


6  *  1,  0  <  r  <  m-1.  Hence,  by  induction  (-l)k£  * 

mr  —  —  m 

m  k 

y  (0  )  ll x)  where  ©  =  (6  )  is  the  associated  lower  trian- 

— ■  mr  r  **=  mn 

o 

v 

gular  matrix  and  ©  is  non-negative  for  all  k.  Hence 


iC’iol  <  z  </> 

r=0 


mr  ’ 
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Lemma  2 . 

A 

<r>  i 

(-D 

kl {k)  (0)  < 
m  '  — 

,m+k. 

1  k  ‘ 

Proof: 

CO 

E 

o 

m  a  t 

u  £  (t 
in 

)  = 

(1  -  u)"1 

exP{-S5'T  (1 

+  a)/ 

(1  - 

u)  )  -> 

E  uVkl(o, 

=  ( 

-*>k 

(i  + 

.  k 
u) 

,  k+ 1 

(-*)k  £ 

<*)Ur 

(1  - 

,-k-l 

u) 

m=0 

(i  - 

u)  x 

r=0 

i 

Hence 

(-i)Vk)  (o)=  Y 

m  r«0 

(%> 

k  .k,  .m+k-r . 
'rM  k  ' 

j? 

ii 

k<*)f 

m+r. 
k  ’ 

.  But 

(in)  ^  (■ m^r )  £  (m^k)  ,  0  _<  r  <_  k,  and  Y  (£)  *  2k ,  so  that  the 

o 

lemma  follows. 

Proof  of  Theorem  2.1.  Clearly*  (m+k)  £  (k  +  l)mk,  m  >_  1, 
so  that  |£^k)  (t)  |/mk  <  k  +  1,  and  |^k)  (t)  |/mK  <  k  +  1, 

£  +  K  £rn(x)  + 

1  <  k  <  K,  m  >  1.  Consider  f(-r)  =  T  (f  in  )  — „ —  +  f Tl  (T) 

—  —  —  rr  m  a  o  o 

1  m 


and 


f(T+h)  -  f ( T ) 


-  E 

m=l 


t  K 
( f  m 
m 


V 


x+h 


V (y) dy 
m 


ni 


K 


+  f 


ll T+h)  -  IA t) 
o  o 


By  dominated  convergence,  we  have  f  '  ( t )  =  7  f+  l A t)  and 

u  in  m 
o 

00  4* 

f '  (x+h)  -  f ' (t)  =  TV  f£'(T+h)  -  £'(t)l  ,  so  that  f’(T)  is 

q  m  L  m  m  J 

continuous  and  goes  to  zero  as  t  -*■  «.  The  argument  may  then  be 
repeated  for  all  derivatives  up  to  and  including  the  K'th.  deriva¬ 
tive.  D 


*  /m+k\  ym+1)  (m-t-2)  ....  (m+k)  <  ^k  2»  3.  .  .  .  (k+1)  _  (^+1)^ 

Iky  l*2»...k  —  1*2. ...k 
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A  simple  extension  of  the  theorem  will  be  of  interest  sub¬ 


sequently. 


ao 

Definition  (a_)  e  S„  <*=>  ^  |  a_  I  <  00 

-  n  K  £  n 

00  ,  . 

Theorem  2.4a.  Let  f(x)  =  Y!  fT  £  (x)  ,  and  let  (f_ 
-  rf  n  n  n 


)  e  S, 


Let  xp  {~r)q  f(T)  =  £  n(T);  P  +  q#<qK.  °  Then 

at  o  n  K  >  1 


’  e  SK-(p+q)  ' 

Proof.  The  theorem  is  proven  by  showing  that  differentiation 
lowers  K  by  at  most  one,  and  multiplication  by  x  lowers  K  by 
at  most  one. 

Thus  let  g(x)  *  £g^n(T),  with  H^lg^l  c  00 »  K  >  1. 

00  00  , 

Then  g(x)  e  C*(R+)  and  h(x)  =  g' (x)  =  [  gj  x)  n<T> 


00 

vhere  | h+ }  £  £  |g  |  ,  as  in  the  proof  of  Lemma  2.2.  Then 


n 1  -  *•*  '"m 
n  m=n 


E"k-X i<  i  i^1  i«ii-  i  i9ii  t  n"-1 1  s  n4i"K< 

o  n=0  m=n  m  m=0  n=0  m 
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Similar  ly  let  s(t)  =  Tg{x)  =  E  g!  (t)  with  (g+)  e  Sv, 

K  >  1.  Then  s(t)  =  £s*^(t)  =  £g*  ((2n  +  l)in(T)  -  (n  +  D^n+1(T)  - 

n^«-i  "•  so  that  (s'!")  is  summable.  Then  s+  *  (2n  +  l)q+  - 

11  1  n  n 

00 

“(n+1)gn+r  ngn-l  '  and  ^  nK_1 1 s^  (  <  »  as  needed. 

o 

Corollary  2.4b.  HnK  |  j  <  J  ~>  TP(§£)q  f(x)  e  C*~P~q  (R+) 

K>1,  P,q^0 
p+q  <  K  ) 

Let  us  now  see  how  smoothness  in  f(t) 
convergence  in  (f*)  .  The  Laguerre  function 
the  second  order  differential  equation 

t  *”(t)  +  r  ( t)  +  (n+Js-*tUn(t)  =  0  (1) 

for  0  <  t  <  «°  .  This  may  be  written  in  the  form 

[4t  -  It  ‘  St]  V*’  ‘  [**  '  D]  £n(t)  (2> 

*  L  V*> 

=  (n+Jj)  *n(t) 


induces  rapidity  of 
fn(t)  satisfies  [  l  ] 
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It  is  easily  seen  that  the  operator  L  is  self-adjoint  in  the 
sense  that 


f  g  (t)  jl,h  (t)J  dt  =  J  ^Lg(t)J  h(t)  dt 


under  simple  conditions,  (4a)  and  (4b)  given  below.  Indeed  one 
may  employ  the  identity 

g  (t)  Dh  (t)  -  h  (t)  Dg  (t)  =  h(t)  Lg(t)  -  g(t)  Lh(t) 

-  *?  [« *  i  -  h  ‘  i?] 

and  integrate  by  parts  to  obtain  (3)  provided  that  the  function 

g  t  ^  -  h  t  vanishes  at  t  »  0  and  <*  and  the  integrands  in 

(3)  are  in  1^(0,°°).  Thus  (3)  is  assured  whenever  h(t)  =  l^it)  , 
n  >  1,  and 


g  1 (t)  is  continuous  and  bounded  on  (0,®):  r  =  0,1,2 


One  further  condition  on  g(t)  is  needed  as  will  be  seen  in  a 
moment.  This  is  the  condition 


Lg (t)  =  fctg(t)  -  tg"(t)  -  g*  (t)  e  L2(0,°°) 


/•°0 

From  (2),  (3),  and  (4a)  we  have  for  g^  ■  I  g(t)  <£  (t)  dt, 
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*  [jf  -  Vt}  dt]  (n+Jj) _1 
**  (n+%)-1  J  ^Lg(t)J  ^n(t)  dt 


00 

If  g(t)  satisfies  (4b),  then  Lg(t)  =  V  (Lg) +  l  (t)  and 

o  n  n 

(Lg)^  *  o(l).  Hence  from  (5)  *  (n+*s)  “:Lo  ( 1)  =  ofn'1). 

The  operator  L  —  (Ht  -  D)k  is  also  se- f-adjoint  under 
similar  conditions  on  the  boundary.  From  (2),  one  has 

Lk  fn(t)  =  { n+  *5 ) k  *n(t)  (6) 

If  both  sides  of  (6)  are  multiplied  by  g(t),  successive  integra¬ 
tions  by  parts  gives 


(n+h)k  9 n  =  /  ^9^  dt 

*u 


provided  that  g(t)  arid  Lmg(t)  satisfy  conditions  (4a, b)  for 

m  **  1»2,..., k-1.  These  conditions  will  be  assured  if  g^r^  (t)  is 

continuous  and  bounded  on  (0,®)  for  0  <  r  <  2k.  We  again  re— 
k 

quire  that  L  g(t)  c  L^fO,®),  and  then  have,  as  above, 

gn  *  o(n  )  ,  n  -*•  (8) 

These  conclusions  are  summarized  in  the  following 
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Theorem  2.5.  If 


a)  f^(t)  is  continuous  and  bounded  on  (0,®),  r  =  0,1,... 2K 

b)  Lrf(t)  e  b2(0,»),  r  =  0,1,. ..,K  ,  (9) 

then 

=  o (n  ^)  ,  n  *►  “  (10) 

n 

It  may  be  noted  that  (9b)  requires  in  general  the  existence  of  a 

2Kth  moment  for  f2(t),  i.e.  that  f  t2K  f2(t)  dt  <  ».  Thus  the 
—  Jo 

coefficients  f^  decay  quickly  to  the  extent  to  which  f(t)  is 

smooth  and  localized  about  t  ®  0. 

Of  particular  interest  is  the  situation  where  f(t)  is  real 

analytic  on  (0,°°)  and  falls  off  quickly  enough  at  ®  so  that  (10) 

can  be  true  for  all  K.  The  class  c”^)  of  "rapidly  decreasing 

functions"  of  interest  to  Fourier  transforms  (see  Dym  &  McKean,  15), 

p.  87)  is  defined  by:  f  (x)  e  C+(R1)  >  x  (^)  f(x)  is 

bounded  on  R1  for  all  non-negative  integers  p,q.  We  may  speak 
correspondingly  of  the  class  c“(R+)  of  rapidly  decreasing  functions 

on  *+-(0,»)  for  which  xq  (§^)P  f (x)  is  bounded  on  R+  for  all 
non- negative  integers  p,q.  Then  as  a  corollary  to  our  previous 
theorem  and  of  the.  earlier  Theorem  2.4  we  have 


-*  0 


Theorem  2.6.  A)  f(t)  e  C+(R+)  if  and  only  if  B)  n  f^ 
n  -*•  «*>  for  all  positive  integers  K. 

Definition 

t  OO 

When  (B)  holds,  (fn)  is  said  to  be  a  "rapidly  decreasing 
sequence"  [5]. 

Note  that  (B)  implies  that  nKAf*  =  nK(f^+^-  -*•  0 ,  n  -*•  <*> 

and  indeed  that  n^Dpf^  ■*•0,  n  where  D  is  any  forward,  cen¬ 

tral  or  backward  difference  operator.  The  name  rapidly  decreasing 

rr  — 

sequence  is  therefore  appropriate.  Such  functions  as  n  0  and 
K  n2 

n  0  are  rapidly  decreasing,  for  0  <  e  <  1. 

It  is  important  to  note  that  smoothness  of  f(t)  alone  does 
K  + 

not  assure  that  n  f^  +  0  for  all  K.  As  an  example,  let 
f  (t)  *  l/(l+t)  e  I«2 .  Then  f  ^  (t)  -*-0,  t  •+■  °°.  We  have 

£1  -  jf  T7t  dt 


e~  9  ( 1+ 1) 


£n(t)  de  dt 


h  r  -1 , ,  -1,  n  -u  . 

-el  u  (1-u  )e  du 
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For  n 

For  1 

Hence 


e  e1*  f  {  (1  -  u)n  u"1  e‘1/u>  du 
*'o 

*  I 

n 

*  2k,  I  =  {f  {  }  du  +  f  {  }  duje”5 

•'o  */l 

>  e*5  f  (1  -  U)2k  u"1  e~1/u  du 
•'l 

1  u  <  2,  £e*  e“1/u  u-1]  >  i  ,  hence  I2k  >  \  J  x2k 

fn  is  not  rapidly  decreasing. 
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7 .  Calculation  of  the  Laguerre  Coefficients  via  Taylor  Series. 

We  have  seen  that,  when  a(t)  e.  L~,  (u)  =  a(i  i^)  is 

regular  about  u  *=  D.  There  are  two  different  Taylor  series  of  value 
in  calculating  the  Laguerre  coefficients.  The  first  valid  in 
principle  for  all  a(i)  in  i2 (R+)  proceeds  as  follows:  We  may 
write 


1+u 

1-u 


) 


+ 


l-u; 


a  ( t)  e 


-t/2 


e-35T(u/l-U)}dT 


The  function  6 (u)  -  u/{l-u)  is  regular  at  u  =  0  and  a(^  +  6) 
is  regular  at  5=0.  Thus 


+ 


u 


(1-u) 


(7.1) 


£ 

k=0 


r  (_a-)k 

Vl-u' 


But 


(— )k 
'l-u' 


Hence 


m 


,m-l,  m 
<k-l)u 


Sc  Ci5 


(7.2) 


(7.3) 
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and  similarly,  since 


u 


(1-u) 


k+1 


<*”> 


(7.4) 


(7.5) 


Thus  one  has 

Theorem  7.1.  If 
are  given  by  (7.3)  and 
Example:  a(t)  = 


a  { t )  r.  12  (H+)  ,  its  Laguerre  coefficients 
(7.5)  . 

t'V6*,  B  >  0  has  ck  =  (^>k  jH  tro+V(e+Js) 

(m+k)  i  <-~)k  (e+Js)  and 


an+  =  I  (-|)k(B+Js)”m“k"1(k)  (m+k)  1 
n  k=0  K 

Note  that  Theorem  7.1  is  useful  when  the  ^  coefficients  are 
available.  In  some  circumstances  a  Taylor  expansion  of  a(s)  about 
8=0  is  more  convenient  because  the  Taylor  coefficients  there  are 
available.  Thus  if 


a(s) 


a (T)e“STdT 
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“a  <5  +  ir^>n  -  „£  V*  . 


0  <  1 u I  <  1. 


k=0 


Thus  from. the  Binomial  expansion  and 


oo 

-  X 


(1-u)  k=  r 


-  r  >  1 


we  have  gnQ  =  (|-)n  and. 


’nk 


2^  (")  (|)n_r  (£l*>  /  k  £  1 


r=l 


_T* 

When  a(x)*efT>0,  A  =  ®,  and  the  coefficients  may  be 
obtained  in  this  way.  Similarly,  when  a(x)  =1,  0  5  t  <  1,  the 

calculation  is  easy. 

Theorem  7.2.  Let  a(s)  be  regular  for  0  £  (s)  <  A,  A  > 


Let  a(s) 


=  ]£ansn,  and  let  ^  +  “(l-u)'1]8  *= 


Then 


angnj.  is  absolutely  convergent. 


n 


proof.  By  Cauchy's  Inequality,  | aR |  <  Mj/R^  where 
•j  <  <  A.  Also 


max  | a (s) | , 

|  s  I  — R 


7-4 


-37- 


By  taking-  6  sufficiently  small,  we  may  assure  that  max  Ik  +  I 

6  1  1-u ' 

<  1*2#  R2  <  and  the  theorem  follows. 
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Difficulty  arises  when  a(t)  and  b(t)  are  known  numerically 
from  data,  say.  Numerical  solution  of  (1)  is  most  naturally 
attempted  by  discretizing  the  continuum,  but  there  are  often  serious 
difficulties,  as  we  will  see. 

The  numerical  solution  of  the  lattice  deconvolution  equation 

#  # 

(2)  is  in  principle  trivial  when  the  coefficients  an  and  bn  are 
known.  For  one  has 


so  that,  recursively, 

f»  .  [b*  -  V  £ V  1  (a*)-1  <4> 

n  L  m=0  m  n_inJ  0 

4 

provided  ^  0.  Now 

00  oo  4-/0 

a?  =  a*  «  /  a(t)  £  (t)  dt  =  /  a(t)e_t/^dt  (5) 

u  u  0  u  0 

ft 

and  one  may  expect,  in  general,  that  aQ  /  0.  One  could  have 
a*  «  0  if  a(t)  is  in  the  space  orthogonal  to  £^(t). 

The  deconvolution  equation  (1)  has  start-up  difficulties  in 
discretization  approaches  when 

a(t)  *  Ktm  /'  t  -*■  0+  ,  K  >  0  (6) 
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for  some  positive  integer  m.  Note,  however,  that  this  does  not 
cause  trouble  in  the  Laguerre  procedure,  as  may  be  seen  best  from 


an  example.  If  a(t)  =  te~at,  a(s) 


(s+a)  ^  and  T* (u) 
21 


(a  + 


1 

1  'l+u.2 

2  I^u) 


|  «  1 

with  a^  =  Ta(0)  =  - -*■  >0.  So  that  there  are  no  start-up 

(a  +  Jj)  " 

difficulties. 

In  those  instances  where  a^  =  0 ,  0  £  m  <_  M-l,  one  will  also 

have  from  T*(u)  =  T*(u)T*(u)  that  b#  =  0,  0  <  m  <  M-l.  Then  T*(u)/T#(u) 

— M  — M  I  ^ 

u  T^(u)/u  Tg(u),  and  the  lattice  deconvolution  proceeds  normally. 

When  a(t)  and  b(t)  are  integral  functions  with  Laplace 
transforms  regular  at  infinity  of  order  at  most  one,  it  may  be  de¬ 
sirable  to  represent  a(t)  and  b(t)  in  terms  of  the  Erlang 
functions  of  Section  4.  This  procedure  has  the  advantage  that  the 
Erlang  building  blocks  are  naturally  ordered  by  their  behavior  at 


t«0 

in  that  tme 

is  smaller 

near  zero 

than  tne-t/^ 

when 

m>n. 

Note 

that  <Ms)  = 

6-(s)/a(s) 

is  regular 

at  ®  when 

a  (s) 

and 

6  (s) 

are  regular 

at  «*>  and 

b  (t)/a(t) 

is  bounded  near  zero. 

The  condition  is  required  for  the  existence  of  any  solution  by 
reasoning  similar  to  that  of  the  previous  paragraph,  since 


f  b  s*(m+1)  /  V  -  _-  (m+1)  _  /u  ,_-(M-N) 

rn 


7  £  a_s 


m 


'  <W* 


-*•  00  where  a  ,  fa¬ 
ro  N 


O  O 

are  the  leading  coefficients  and  M>N  is  required.  The  solution 
f(t)  will  then  itself  be  in  the  Erlang  family,  and  lattice  decon¬ 
volution  of  the  Erlang  coefficients  is  viable. 
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9.  Numerical  Examples  of  the  Kethod 

We  present  some  numerical  examples  of  the  transform  methods 
for  the  calculation  of  multiple  convolutions  as  well  as  for  the 
solutions  of  Volterra  equations  (0.1),  (0.2).  All  computations 
were  done  on  a  Burroughs  6750  computer  in  a  time-sharing  mode  using 
APL  as  a  programming  language.  The  computations  were  in  single 
precision  (11  digits) .  Relevant  formulas  were  usually  coded  in 
an  "obvious",  straightforward  way,  with  no  attempt  made  to  opti¬ 
mize  the  subroutines  for  speed  or  accuracy.  In  spite  of  this,  the 
results  displayed  in  this  section  were  typically  gotten  with  CPU 
times  measured  in  seconds,  and  with  no  evidence  of  numerical  prob¬ 
lems  such  as  roundoff,  truncation,  underflow,  etc. 

The  Laguerre  functions  were  calculated  using  the  recursion 
relation  (see  Appendix  A) 


- f2  - 


(t+1) 

n+1 


un(t)  -  [ 


n 


n+lUn-l (t) 


with  £Q(t)  =  exp(-vt)  and  £^(t)  =  ( 1— t )  exp  (-vt)  .  Similarly, 
e  (t)  =  (^)e  i (t) ,  with  eA(t)  -  exp(-vt).  The  use  of  the  scale 
factor  v  is  recommended  for  computational  flexibility.  When  v  =  %, 
the  above-mentioned  Laguerre  and  Erlang  functions  are  obtained.  When 
v  *  0,  the  Laguerre  functions  become  the  Laguerre  polynomials  and 
the  Fourier-Erlang  series  reduces  to  the  Taylor  series.  The  use  of 
the  scale  factor  v  'is  equivalent  to  the  use  of  the  exponential 
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then 


transform  of  Section  5,  for  if  fQ(t)  =  f(t)e 

£„(t)  =  I  f+e~vtL  (t)  for  v  =  e  +  h.  Hence  the  use  of  v  *  'j  is 
e  n  n 

justified  providing  f(t)  is  Laplace  transformable  and  f 0 (t)  e  L2- 
Typically,  the  computer  will  give  the  same  answer  for  any  v  in 
the  vicinity  of  h,  and  the  user  may  choose  a  value  of  v  empirically 
to  get  the  best  numerical  accuracy. 

The  Laguerre  and  Erlang  coefficients  for  exp (-at)  are  needed 
frequently  in  applications  and  are  derived  here  for  easy  reference. 
For  f  (t)exp  [  (v-*s)  t]  £  L2(0,«),  (1.5)  becomes 


f^+(v)  =  /  f  (t)  exp{  (v-!j)  t}fn  (t)dt 


(1 


(in  which  £n(t)  is  the  classical  Laguerre  function)  and  for 


f(t)  =  exp (-at),  we  recognize  the  Laplace  transform  of  £n(t)  evalu¬ 


ated  at  a  +  h  -  v.  From  (1.8), 


f“<v)  -  , 


n  =  0,1,2, . . . 


(2 


(a+l-v) ' 


For  the  Erlang  transform  we  have 


f(t)  =  exp(-at)  =  l f^+ (v)  (tn/n! )exp(-vt)  , 


so  that  l  f^f (v) (tn/n! )  =  exp(v-a)t,  and 


f?+(v)  *  (v-a)n 


n  =  0,1,2, .. . 


(2 


-<  3- 


These  coefficients  for  the  exponential  function  serve  as  useful 
building  blocks  for  calculating  coefficients  of  more  involved  func¬ 
tions.  For  example,  if  f(t)  is  an  Erlang  density,  its  coefficients 
may  be  easily  computed  by  recognizing  f (t)  as  an  integral  number  of 
convolutions  of  an  exponential  density. 

Example  1;  Server  Busy  Period  Density 

This  example  illustrates  the  use  of  the  transform  methodology 
for  the  evaluation  of  multiple  convolutions .  Reference  [3]  illus¬ 
trates  the  numerical  effort  that  may  be  required  to  perform  multiple 
convolutions  by  numerical  methods.  When  the  function  being  convolved 
has  rapidly  decreasing  coefficients,  our  method  yields  a  fast  and 
accurate  algorithm. 

Consider  a  single  server  M/G/l  queueing  system  [18]  having  a 
Poisson  stream  of  customers  at  rate  x  and  service  time  distribution 
which  is  absolutely  continuous  with  density  a(t).  The  server  busy 
period  density  s(t)  is  well  known  [18]  and  is  given  by 

s  ( t )  =  l  exp(-Xt)  |rc+lTi  a^n  +  1^(t)  *  (3) 

The  study  of  s(t)  has  been  hampered  by  its  relative  numerical  intract¬ 
ability.  Only  when  a(t)  is  a  pure  exponential  is  a  closed  form  answer 
available.  We  take  a(t)  to  be 

a (t)  -  2  Ce_t  -  e'2t)  (4) 
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i.e.,  a{t)  is  the  density  of  the  sum  of  two  independent  exponentially 
distributed  variables,  one  with  parameter  1,  the  other  with  parameter  2 
The  Laguerre  or  Erlang  coefficients  for  a(t)  are  written  down  imme¬ 
diately  using  (2a)  or  (2b) .  By  inspection,  we  see  that  the  best 
value  of  v  for  the  Erlang  coefficients  is  1.5,  and  that  the  best 
value  of  v  for  the  Laguerre  coefficients  will  be  between  1  and  2 
(turned  out  to  be  about  1.3). 

To  perform  multiple  convolution,  we  make  use  of  (1.13)  and 

(4.6),  which  state  that  for  either  Laguerre  or  Erlang  building  blocks 

we  have  T*  *  T*T#.  Hence  tJ  =  (T?)m,  where  h(t)  «=  a {m)  (t)  and  we 
i  y  X  g  n  8l 

have  the  following  algorithm  for  the  multiple  convolution  of  a (t) 
using  Laguerre  building  blocks: 

I.  Generate  or  store  in  the  computer  the  Laguerre  coefficients 
+  N 

{a  )_  for  a(t).  Note  that  this  is  a  finite  set. 
n  u 

II.  Convert  {a+}  to  {a*}.  Since  (u)  =  (l-u)T*,  this  corre- 
n  n  a  a 

sponds  to  a  simple  differencing  operation. 

^  N 

III.  Perform  m-fold  discrete  convolution  on  (aw}„  .  Retain 

n  o 

only  the  first  N+l  terms  in  each  convolution.  The  result 

#  N 

is  the  sequence  { h^ J  q . 

s  4* 

IV.  Convert  {fr}  to  {h  }.  This  is  the  inverse  of  the  differ- 
n  n 

encing  operation,  i.e.,  summation. 

N 

V.  Sum  the  series  I  h  (  (t)  to  get  h(t). 

0  n  n 

A  completely  analogous  procedure  is  used  for  the  Erlang  building  blocks. 
A  few  remarks  may  be  helpful: 
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(1)  The  steps  involve  a  finite  number  of  computer  operations 
on  discrete  sequences.  The  operations  are  things  that  a  computer 
does  well.  The  steps  are  easy  to  program  and  are  typically  done 
with  high  speed  and  accuracy. 

(2)  In  applications,  step  III  was  done  using  the  formula 

*  |  Vn-k  ' 

Since  N  is  typically  less  than  100  and  frequently  less  than  50,  the 
use  of  sophisticated  techniques,  e.g..  Fast  Fourier  Transform,  for 
the  discrete  convolution  does  not  seem  warranted. 

(3)  A  typical  sequence  of  operations  is  to  make  a  trial  run 
with  N  *  30,  repeat  for  various  values  of  v  to  get  the  most  rapid 
convergence,  then  increase  N  as  needed  to  get  the  desired  accuracy  at 
all  necessary  values  of  t.  Thus,  working  with  the  computer  in  an 
interactive  mode  is  ideal  for  this  methodology. 

With  this  algorithm,  the  functions  a^n+^(t)  may  be  computed 
for  desired  values  of  t,  starting  with  (t)  =  a(t).  Figure  1  dis¬ 

plays  s(t)  for  X  =  .2  and  A  =  1.  For  both  the  Laguerre  and  Erlang 
building  blocks,  accuracy  was  best  near  the  origin  and  decreased 
with  increasing  values  of  t.  Three  digit  accuracy,  i.e.,  sufficient 
to  plot  figure  1,  was  gotten  by  using  N  =  30  for  all  the  convolutions 
(either  set  of  building  blocks)  and  using  the  first  6  terms  of  the 
infinite  sum  for  s(t).  When  high  accuracy  was  desired,  the  Erlang 
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building  blocks  proved  superior  to  the  Laguerre  building  blocks. 

Using  N  «  40  for  all  convolutions,  accuracy  to  about  10  digits  after 
the  decimal  point  was  attained  out  to  t  =  10  using  23  terms  in  the 
infinite  sum  expression  for  s(t). 

The  service  time  distribution  has  mean  2/3,  so  that  for  X  =  .2, 
s(t)  represents  an  honest  density,  i.e.,  it  integrates  to  unity. 

For  X  >  2/3,  the  formula  for  s(t)  remains  valid  but  integrates  to 
less  than  unity  since  there  is  positive  probability  that  the  queue 
will  never  empty  after  the  first  customer  arrives. 

Example  2:  Solving  a  renewal  equation 

This  example  illustrates  the  transform  methodology  in  the  solu¬ 
tion  of  a  Vol terra  integral  equation  of  type  2,  e.g.,  a  renewal  equa¬ 
tion.  The  use  of  series  methods  in  the  solution  of  renewal  equations 
has  been  used  with  some  success  by  previous  authors,  e.g.,  [8],  [16], 
who  have  typically  restricted  themselves  to  the  Taylor  series.  As 
this  paper  makes  clear,  the  use  of  the  Taylor  series  is  a  special 
case  of  a  much  more  general  approach. 

In  a  reliability  context,  suppose  an  item  has  failure  time 
density  a(t),  and  suppose  that  the  maintenance  policy  is  to  replace 
the  item  at  random  times,  whether  it  has  failed  or  not.  Let  the 
replacement  times  be  exponentially  distributed  with  rate  parameter 
X,  and  let  s(t)  be  the  density  of  the  resulting  effective  failure 
time  distribution.  Then  probabilistic  reasoning  gives 
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s  (t)  =  b  (t)  *{  6  (t)  +  c(t)  +  c^(t)  +  ...} 

=  b(t)*J  c(k)  (t) 

0 


(5) 


where 

b(t)  =  a(t)exp(-Xt) 

(6a) 

00 

c (t)  =  Xexp(-Xt)  /  a(y)dy 

(6b) 

t 


and  S (t)  is  the  Dirac  delta  function.  It  is  clear  that  s(t)  obeys  a 
renewal  equation 

s  ( t )  =  b(t)  +  c  (t)  *s  (t)  (7) 

The  multiple  convolutions  in  (5)  could  be  evaluated  as  in  the  preceding 
example,  but  the  direct  solution  of  the  renewal  equation  (7)  is  much 
more  economical  and  accurate. 

To  solve  (7),  we  note  that  (7),  (1.12),  and  (4.5b)  imply  that 


«  #  #  # 

T  =  T  +  T  T* 
s  b  c  s 


which  implies  the  discrete  renewal  equation 


(s?)  =  {b?}  +  {c?}*{s?} 
n  n  n  n 


(8) 


i.e. , 


s*  «  bf  +  l  c^s'  , 
n  n  g  k  n-k 


n 


(9) 


#  |  i 

Equation  (9)  can  be  solved  for  sn  recursively,  starting  with  =  b", 

I 

providing  that  c.  f  1.  The  algorithm  for  solving  (7)  is  hence  steps 
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I  -  V  of  Example  1,  except  that  step  III  requires  the  solution  of 

the  discrete  renewal  equation  (8). 

-t  -  2 1 

We  choose  a(t)  -  2 (e  -  e  )  as  in  Example  1  so  that  (6a), 

(6b)  become 


b(t) 
c  (t) 


2  [e 


2  Ae 


— (1+A ) t 
-  ( 1+  A )  t 


e-(2+A)tj 


Ae 


- (2+A) t 


(10a) 

(10b) 


The  Laguerre  or  Erlang  coefficients  for  b(t)  and  c(t)  may  be  written 
down  immediately  by  using  (2a) ,  (2b) ,  respectively.  The  resulting 
s(t)  is  displayed  in  figure  2  for  A  »=  .  2  and  A  =  4.  As  expected, 
the  mean  time  to  complete  the  message  transmission  increases  as  the 
mean  interruption  rate  increases.  Using  N  =  4-0,  both  sets  of  building 
blocks  give  about  4  digits  of  accuracy  after  the  decimal  point  in  the 
time  interval  [0,  10]  for  A  =  4.  For  A  =  .2  and  N  =  40,  we  get  5 
and  10  digit  accuracy,  respectively,  for  the  Laguerre  and  Erlang 
building  blocks. 

Example  3:  The  deconvolution  problem 

The  deconvolution  problem  involves  solving  an  equation  of  form 

b  (t )  -  a  (t)  *f  (t)  (11) 

with  a(t),  b(t)  known  and  f (t)  unknown.  This  problem  is  discussed  in 
section  8,  and  equation  (8.4)  provides  an  algorithm  for  solving  the 
discrete  deconvolution  equation  (8.2),  which  is  of  the  form 

{bn}  =  {aJ)MfJ}  *  (12) 
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It  should  be  noted  that 
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a.e. 


(12)  can  be  recast  as  a  renewal  equation. 


it,*) 


+  lfi(n)  - 


n  n 


where 


6(n)  =  < 


1  n  =  0 

0  n  ?  0 


Hence,  a  subroutine  for  solving  a  discrete  renewal  equation  also 
serves  to  solve  the  discrete  deconvolution  equation. 

In  an  M/M/1  queue  setting,  let  sn(t)  be  the  density  of  the  time 
for  the  queue  to  empty,  given  a  start  with  n  customers.  It  is  known 
(9]  that  sn(t)  satisfies 


9-n(t)  “  sn(t)*g0(t) 


(13) 


where  gn(t)  -  (A/w)  U+w)  tIn  (2/Xut) ,  n=  0,±1,...  (14) 

with  A  the  arrival  rate  and  y  the  service  rate.  The  density  sn(t)  is 
known  19]  to  be  equal  to  (n/t) g_n (t) ,  but  we  solve  (13)  numerically 
to  illustrate  the  deconvolution  technique. 

When  >  =  y  *=  h,  (14)  simplifies  to 

g  (t)  -  e_tl  (t)  . 

n  n 

We  expand  I^tt)  in  a  Taylor’s  series  and  immediately  get  a  natural 
expression  of  gn(t)  as  an  Erlang  series  with  v  =  1.  The  deconvolution 
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algorithm  proceeds  with  no  difficulty  and  figure  3  displays  s^ft). 

The  density  is  integrable,  but  has  a  very  heavy  tail,  i.e.,  it  has 
no  moment.  Twenty  terms  in  the  series  provides  more  than  enough 
accuracy  to  make  the  plot,  but  more  terms  are  needed  for  large 
values  of  t.  Sixty  terms  gave  about  5  digits  of  accuracy  after 
the  decimal  point  at  t  =  50,  with  10  digits  at  t  *  20. 

When  A  =  .25  and  p  =  1,  (14)  becomes 

g_  (t)  =  (Js)ne“1,25tIn<t) 
n  n 

Again,  the  expansion  of  1^  (t)  in  a  Taylor's  series  provides  a  natural 
expression  of  9n(t)  as  an  Erlang  series,  this  time  with  v  =  1.25. 
Figure  3  also  displays  (t)  for  this  case.  The  tail  is  very  light 
and  30  terms  in  the  series  gives  an  accuracy  of  about  8  digits  out  to 
t  *  10. 
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Fig.  3:  s.(t)  versus  time 


10.  Interpolation  Methods  and  Problems 

For  the  solution  of  integral  equations  of  form  (0.1)  or  (0.2) 
where  the  given  functions  a(x)  and  b(x)  are  only  known  numerically 
at  some  finite  set  of  points,  interpolation  procedures  are  needed. 
Interpolation  of  the  data  via  spline  functions  may  or  may  not  be 
suitable  for  the  methods  we  have  described.  The  discontinuities 
in  the  derivatives  at  the  knots  intrinsic  to  those  methods  preclude 
the  Erlang  transform  methods  of  Section  4  since  a(x)  and  b(x)  cannot 
then  be  integral  functions.  Laguerre  transforms  are  possible,  but 
the  functions  will  not  be  rapidly  decreasing  and  hence,  as  shown 
in  Section  6,  the  Laguerre  transform  coefficients  are  not  rapidly 
decreasing  and  the  tails  of  a^  could  be  excessively  long.  A  suffi¬ 
ciency  of  data  points  could  permit  spline  interpolation  of  high 
order,  and  corresponding  accuracy.  The  difficulties  surrounding 
such  finite  data  set  problems,  and  the  accuracy  of  our  methods  when 
applied  to  such  problems,  require  further  study. 
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Generalization  of  the  Method  to  Other  Families  of  Function 

The  basic  method  can  be  extended  to  other  classes  of  functions. 

Consider  any  sequence  of  functions  with  Laplace  transforms 

(a  <»))•  for  which  a  (s)  has  the  form 
no  n 

a  (s)  =  a  (s)  (6(s))n  ,  (1) 

n  o 

Suppose  further  that  under  the  change  of  variable  w  =  3(s), 

a  (s)  becomes  p(w),  i.e.,  a  (&-1(w))  =  p(w),  where  p(w)  is 

o  o 

regular  at  w  =  0.  Let  g(t)  =  £gnan(t),  h(t)  =  £hnan(t).  Then 

l[g(t)]  *  £gn°n(s)  =  ao{s)  Hgn[eo(s)]  n  =  p(w)  and  L[h(t)]  * 

p(w)  J  hnwn.  In  our  previous  notation,  we  then  have  for  T^(w)  =  £gnwn 

p(v)  Tg*h  -  p(w)  Tg  (w)  p(w)  T*(w).  (2) 

Hence 

Tg*h(w)  "  Tg(w)  Th(w)  (3> 

#  t 

where  T  (w)  *  p(w)  T  (w) and  the  framework  for  computation  developed 
g  g 

in  Section  1  for  the  Laguerre  transform  carries  over. 

There  are  many  classes  of  functions  ia^Ct)}  having  the 
necessary  algebraic  structure  for  a^fs).  Two  such  classes  quite 
different  in  form  from  the  Erlang  class  and  Laguerre  class  are  of 
interest.  In  the  first  the  building  blocks  are  the  functions 
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an(t) 


4  TT  t 3 


exp  {-nV4t} 


which  arise  in  the  theory  of  diffusion.  For  these  functions  one 


-n/s~ 


an(s)  =  e 


B(s)  =  e  ,  ao(s)  =  = 


The  second  family  has  Bessel  function  building  blocks 


V*’  *  Jn(t>  ' 


for  which  II] 


an(s)  *  {  Vs2+1  -  s}n/  f42+l 


b (s)  =  { ys*+i  -  s>  (7b> 

aQ(s)  =  (s2+l)  **  ;  .  p(w)  =  2w/(l+w2).  (7c) 

For  any  of  these  families  (a  (t))  we  may  use  instead  (a  (t)e-^T), 

n  n  ' 

since  the  necessary  properties  are  preserved. 

In  special  contexts,  such  as  underwater  sound  signals,  a 
particular  family  might  be  more  appropriate  than  the  Erlang  family 
or  Laguerre  family.  Each  family  of  functions  has  special  advantages 
and  special  limitations  and  computational  algorithms  appropriate 
to  that  family  are  needed. 
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Inversion  of  Laplace  Transforms 


The  operational  framework  we  have  developed  for  convolutions 
also  provides  formulae  for  inverting  Laplace  transforms  in  terms  of 
the  building  blocks.  Two  such  formulae  of  general  interest  are  pre¬ 
sented  in  this  section.  Other  formulae  may  be  obtained  from  Section 
11  for  special  classes  of  functions. 

OD 

When  a(s)  =  /  e~sta(t)di  is  regular  at  infinity  and  vanishes 
0 

there,  inversion  of  a(s)  as  a  series  is  implicit  in  Theorem  4.1. 

This  may  be  worth  stating  formally.  The  result  is  given  in  a  more 
restrictive  setting  by  D.  V.  Widder  [20] . 

Theorem  4 ♦ 3 

If:  (a)  a (s)  is  regular  at  ®  and  vanishes  there 

(b)  =  i  <=m(;)m+1 

then  a(-r)  is  entire  and  given  by 


a  (  t  )  = 


I 

0 


m 


m 

T 

m2 


>-t/2 


where  the  series  is  absolutely  convergent  for  all  t. 

Note  that  the  Taylor  expansion  of  TL  (u)  about  u  =  0  together 
with  (1.4)  also  provides  an  explicit  inversion  formula  for  Laplace 
transforms.  Specifically,  one  has: 

Theorem  4.4 


Let  a  (  t  )  c  L«(0,»)  and  let  a(s)  •»  i[a(r)].  Then  a(t)  is  given  by 
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(13) 


a(t) 


*  I 

m=0 


,1  ,d  smr  1  ,1  l+u>, 

{5T(d[I)  HI*1 


u=Q 


£m(T)  } 

m 


Numerical  inversion  of  Laplace  transforms  via  numerical  inte¬ 
gration  techniques  has  received  considerable  attention,  e.g.,  in 
the  book  by  Bellman,  Kalaba,  and  Lockett  [2].+  Such  techniques  do 
not  appear  to  exploit  analytically  founded  accuracy  as  fully  as 
our  procedures  and  are  placed  in  jeopardy  by  ill-conditioning  [2, 
p.  33] .  An  extensive  comparison  of  the  two  procedures  is  needed. 

The  idea  of  inverting  Laplace  transforms  via  the  operational 
tools  of  Laguerre  functions  goes  back,  to  the  authors'  knowledge, 
to  Lanczos  in  his  book.  Applied  Analysis  [11] .  The  presentation 
there  is  somewhat  obscure  and  the  scope  of  the  procedure  and  its 
legitimacy  are  essentially  undeveloped.  Feller  presents  explicit 
inversion  formulae  for  Laplace  transforms  [6]  of  largely  theoretical 
value.  Neither  Feller  nor  Widder  nor  Bellman  refer  to  Lanczos' 
method. 
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+The  reader  may  wish  to  examine  the  papers  by  A.  Papoulis,  "A  new 
method  of  inversion  of  the  Laplace  transform,"  Q.  Ap.  Math.  14, 

1957  and  R.  V.  Churchill,  "The  inversion  of  the- Laplace- transforma¬ 
tion  by  a  direct  expansion  in  series,"  Math.  Zeitschrift  42,  1937. 
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Appendix  A 


Classical  properties  of  the  Laguerre  polynomials  L  (x) 


and  functions  £n(x)  =•  Ln(x)e-X//2  (cf.  f  1  ] ) 


n 


1.  /  e  XL  (x)L  (x)dx  = 

n  m 


=  6 


mn 


2.  J  *n<x)  £„,<*> dx  - 


3.  Y  unLn(x)  =  (l-u)“1exp{- 


A. 


I  un£n(x)  =  (l-u)”1exp{-i  x 


n 


5.  L(x)  =  l  (-l)k  iy  (J)xk 
n  k=0 


n 


6.  xn  =  I  (-l)k  £y  (j?)Lk(x) 

k_0 


7. 


/  e’^UJdx  -  5^  (|^)n  ;  /  2nU)dx 


(-1) n • 2 


8.  L  (0)  «  £  (0)  =  1 

n  n 


9.  £  (x)  ^  e"x/2  ,  x  -  » 

n  n  I 


t  1  x  id  >  n  /  n  "X  i 

10-  Vx>  ■  ft  e  'a?1  lx  e  > 


«•  I  Vx)  ST  *  eue_1*xJA(2^xu) 
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12.  xln(x)  =  (2n+l)£n(x)  -  (n+1) *n+1 (x)  -  n^n_1  (x) 


13.  -£n(x)  =  h£n(x)  +  I  *_(x> 


A  TABLE  OF  LaGUERRE  TRANSFORMS 


THIS  PAGE  IS  BEG  I  Q 

VHuil  CGi'I 
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explanat;;-'  of  some  entries  in  the  table  of  operational  properties  of  laguerre  transforms 

(Note] 

Tot  notational  convenience,  any  element  with  e  negative  subscript  is  considered  to  be  aero 
in  the  following  discussion. 


♦  l  1^  *  T  X!* 

oxt).  \  -js  /  a"  '  1  ~  -fr- 

/  \  a* 


Lr— F; 


fe**,  *  £,f:  -  -3L  4?  -  XlK-X*:]  R-CO 

«tw  •  s*w  • 

of  -  a!  ~  ^  £  -  21  ■ $1  =  "T  £*  -  f  W 

ft?-  a:  -  dl.,  /!*$;_  f,C  -t-X  C  -  t  >  V-  •‘*'1 


aCc>-T$«)  v.  /  al  “  “  Cn*°  ‘ 

^)feL,lo.0»]  /  \  q*  *  -  CV"-»?£L  TCv»-oCa 

•o  ••  ^ 

act)  »  T^-tt)  *  X^  X4-.AX)  *  XC  -  C«r.)  S^OtJ  _  «  jL„fctA 

MX •  MS* 

«  X  -  c«ti)  ^2c.  -  *  C-i  1  *-~Ct) 

M«£> 

*•<  «C-£ 

oZ  ~  (£  -aI-,  *  U*"-roC  -  C«"tv) 

*  -  csw- «>  fl,  ■*■  c.*-o  f*.v 


atx)-^  foi^o  \  / 

£.if;i<  o°  /  ^  \  «.*.  -  — sf*  t  ->  £>o-C_ 


(r**F ) 


dtt)  »  [  -fwA)  -  .-.  «<0O  =  X  -  <»w] 

T.»f>  -  -  {-&&)  -  -^=f  bM  -  *  kS&5>]  -  ^  tV<-">  -  Tf«»] 

■C’1"'  =~rhr  t-C  -z  ii'l^rv 

••• «:  -  -*  £<->'■”  C  -  C~ 

*»«i 

dl  .  d  -d.  -  -i  £eo-[  e„.  -  C-l  -  -t .[T*-i  Z>’“  C-] 

-  -  i  f*  <■  ^  5>>-  C. 
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•  v  n  *  s  n J 


0.  Introduction  and  Summary 

In  a  previous  paper,  hereafter  designated  by  [A],  a  description  was 

given  of  a  Laguerre  transformation  which  maps  a  function  f(r)  in  L2(0,“) 

$  00 

into  a  sequence  (f^Q  on  the  nonnegative  integers.  Moreover,  for  two 
such  functions  f(T),  g (t ) ,  the  convolution  f(x)*g(t)  is  mapped  into  the 

#  r  S  #  #  # 

lattice  convolution  (f*g)  =  l  f  g  =  ( (f  )*(gm))  •  One  obtains  thereby 

n  a  ii— m  in  m  m  ri 

u 

an  algorithmic  basis  for  the  computation  of  multiple  convolutions  f  (x) 
and  related  infinite  series  of  importance  to  statistics  and  applied  proba¬ 
bility. 

Such  Laguerre  transforms  have  one-sided  functions  as  their  natural 

domain  because  the  Laguerre  polynomials  (t)  and  Laguerre  functions 
“  T  /  2 

£n(r)  =  Ln(x)e  are  associated  with  the  one-sided  weight  function 
e~T  on  (0,°°).  Nevertheless,  the  methods  have  a  simple  extension  to  two- 
sided  functions  on  the  full  continuum  (-*,  <*>)  via  the  same  Laguerre 
functions  as  we  will  see. 

A  variety  of  applications  exist  to  statistics,  operations  research, 
and  engineering.  In  statistics,  for  example,  one  has  need  for  multiple 
convolutions  of  two-sided  distributions  unavailable  analytically,  that 
of  the  logistic  distribution,  for  example.  Even  relatively  innocuous 
distributions  such  as  the  Laplace  distribution  convolve  with  difficulty. 

In  operations  research  studies  dealing  with  queues,  inventories  and 
storage  systems,  one  encounters  as  a  structural  entity  [3]  the  extended 

oo 

r  00 

renewal  density  h(x)  =  )  a  (x) ,  where  a(x)  is  a  probability  density 

1 

function  with  two-sided  support.  For  many  densities  of  interest,  evalua¬ 


tion  of  h(x)  has  been  resistant. 
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In  the  earlier  paper  [A],  the  crucial  role  of  the  complex  plane  in 
the  formulation  of  the  algorithms  was  evident,  even  though  the  algorithms 
were  entirely  in  the  real  domain.  For  the  bilateral  transform,  the  com¬ 
plex  plane  is  again  very  much  present,  with  Laurent  expansions,  bilateral 
Laplace  transformation  and  conformal  mapping  entering  as  crucial  tools. 

The  first  section  extends  the  earlier  formalism  to  the  full  continuum. 
That  this  extension  is  natural,  and  not  just  an  artificial  piecing  together 
of  the  formalism  for  each  half-line,  will  be  clear  from  (1.9),  (1.12)  and 
(1.13).  The  harmony  of  the  basis  will  also  emerge  vividly  in  Section  3, 
which  deals  with  the  extent  of  the  transform  coefficients,  and  associated 
uncertainty  relations.  The  topic  of  extent  is  crucial  to  the  utility  of 
the  Laguerre  transform  method  as  a  numerical  tool.  Numerical  examples  are 
presented  in  Section  5.  A  table  of  contents  provides  the  reader  with  an 
overview  of  the  paper. + 


+Two  references  (V.  I.  Krylov  and  N.  S.  Skoblya  [8],  and  W.  T.  Weeks  [12]) 
have  come  to  the  authors'  attention  subsequent  to  publication  of  [A].  Both 
deal  with  the  use  of  Laguerre  functions  for  the  numerical  inversion  of  one¬ 
sided  Laplace  transforms. 


1.  The  bilateral  Laguerre  transform 


In  this  section  the  one-sided  methodology  of  [A]  will  be  recast  in 
a  natural  way  to  provide  a  corresponding  representation  of  two-sided  func¬ 
tions.  In  this  extended  setting  the  one-sided  functions  will  be  a  special 
subcase. 

Let  f(x)  be  any  function  in  L 2 Let  U(x)  *  1,  x  £  0,  U(x)  =  0, 
r  <  0.  Then  f(x)  =  f+(x)  +  f  (x)  where  f+(x)  =  f(T)U(t)  is  in  and 

f  (x)  =  f(x)U(-x)  is  in  L2(-°°,  0).  The  discrepancy  at  x  =  0  may  be 
ignored.  Clearly,  one  may  write,  employing  the  notation  in  [A], 

00 

(1.1)  f+(T)  =  l  ff  t  (T)U(T) 

0 

00 

(1.2)  f  (T)  -  l  f l  l  (-T)U(-T) 

0 

and 

00 

(1.3)  f(T)  =  I  f>m(T)  ,  X  X  0  , 

.00 

where 

(1.4a)  hm(T)  =  ^m(T)U(x)  ,  m  >  0 

(1.4b)  hjn(x)  =  -£_m_1(-T)U(-x)  ,  m  <  0  . 

The  puzzling  minus  sign  in  (1.4b)  plays  an  important  role  which  will 


emerge  soon.  We  see  that 
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Note  also  from  (1.1),  (1.2)  and  (1.3)  that 


(l.M)  <  -  C  ,  n  -  0  »  <  -  ,  »  <  0  . 

The  set  of  functions  (h  (t)}_co  form  a  complete  orthonormal  system  for 
l~2  (-•,<*>)  .  Consequently,  one  has 


oo 

(1.5)  f+  =  /  f (t)h  (r)dr 
m  J  m 


/  f(x)£  (r)dt  ,  m  2  0 
0  m 


-/  f  (-t)f.  .  (t)dT  ,  m  <  0 

l  0 


From  the  orthonormality,  one  has  the  Parseval  relation 

(1.6)  J  f2(r)dx  =  l  f*2  . 

-00  -00 

We  know  that  Bilateral -Laplace  transformation  gives  for  n  £  0  (cf.  [A]), 

00 

a. 7)  WT)U(T)1 =  /  e"ST  VT)u(T)dT 

—  00 

=  (s  -  l/2)n/(s  +  l/2)n+1  ,  Re(s)  >  -  1/2 

and 

(1.8)  l-B[-^n(-T)U(-T)]  =  (s  ♦  1/2 )n/  (s  -  l/2)n+1  ,  Re(s)  <  1/2  . 

It  then  follows  from  (1.4a)  and  (1.4b)  that,  for  all  n  (explaining  the  minus 
sign  remark  below  (1.4b)), 

(1.9a)  LB[hn(t)]  =  (s  -  1/2 )n/ (s  *  l/2)n+1  ; 

<  n  <  a,  .  _  1/2  <  Re(s)  <  1/2  . 
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From  (1.4)  and  (1.9)  we  have  the  interesting  identity 


(1.9b)  hm(T)*hn(T)  =  h0('O*hm+n(T)  1  all  m.n.r 

which  may  be  compared  with  eq.  (1.9)  of  [A].  For  f(t)  e  with 

(f^)  €  £j,  i.e.,  'l  | f n |  <  <*>,  we  have  formally  from  (1.3)  and  (1.9a) 


(1.10) 


♦B(S)  d=f  Igtf(T)] 


f+rs  -  1/2, n  1 
Vs  +  l/2}  s  +  1/2  * 


A  discussion  of  the  domain  of  validity  of  (1.10)  will  be  given  in  the 
next  section.  The  ideas  and  results  are  a  two-sided  extension  of  the  corre¬ 
sponding  one-sided  results  presented  in  [A].  As  we  will  see,  the  series  in 
(1.10)  will  be  absolutely  convergent  in  some  strip  -A  <  Re(s)  <  A+,  A_  >  0, 
A+  >  0,  when  f+(x)  and  f_(t)  are  sufficiently  smooth  and  "rapidly  decreasing" 
as  in  [A].  It  then  follows  from  (1.10),  that  in  the  corresponding  domain  D 
in  the  complex  u-plane,  containing  (u:  |u|  =  1),  one  will  have  for  the  dagger 
generating  function  T^(u) 


(1.11) 


Tj(u) 


d«£I 


rt  m  1 
f -u  =  — 
m  1-u 


k  rl  IlHn 
B  2  1-u' 


and 


\ 

* 

I 


(1.12)  Tj(u)  ¥  (l-u)T+f(u)  =  g) 


For  bilateral  convolution  it  is  known  that  Lg[f  (r)*g(t)]  =  4>g(s)ifg(s)  .  Hence 
we  have  from  (1.12), 

(1.13)  Tf*g(u)  =  T*(u)T*(u)  . 

It  then  follows  that  as  in  [A],  bilateral  convolution  may  be  mapped  into 
lattice  convolution  via  the  two-sided  Laguerre  transform 
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f *  =  J  f(T)[hn(?)  -  hnl(T)]dT  , 

.V 

i.e. , 

#  +  + 

(1.14)  f  =  fT  -  fT  ,  . 

v  J  n  n  n-1 

The  validity  of  (1.10)  in  a  convergence  strip  containing  the  imaginary 
axis  in  its  interior  will  permit  T^(u)  to  be  regular  in  an  annulus  contain¬ 
ing  the  unit  circle  in  its  interior,  when  is  regular  at  infinity  as 

we  will  see  in  Section  2.  The  situation  is  somewhat  modified  in  the  absence 
of  such  regularity. 

Note  that  f(-r)  even,  i.e.,  f(t)  =  f(-t)  implies  that  is  even  in 

s.  It  then  follows  from  (1.12)  that 

(1.15a)  f(t)  =  f(-T)  <=>  f*  =  f*n  , 

i.e.,  f(t)  even  in  t  implies  that  f^  is  even  in  n.  Similarly,  one  has 

(1.15b)  f(t)  =  -f(-T)  <=>  f*  =  -f*n  . 

The  reader  will  verify  from  (1.1)  and  (1.2)  that  when  f(x)  is  even 
(so  that  f+(t)  =  f_(-t)),  then  f^+  =  f^  .  It  then  follows  from  (1.4d) 
that 


(1.16a)  f(t)  =  f(-x)  <=>  f*  =  -fj_n_n 

(1.16b)  f(t)  *  -f(-t)  <=>  =  f+ 

n  (-n-1) 

i.e.,  the  daggers  have  their  symmetry  about  n  =  -  1/2. 
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It  is  clear  from  the  simplicity  of  the  form  of  (1.13)  and  from  the 

relative  simplicity  of  the  symmetry  relations  (1.15a,b)  for  f  over  those 
f  # 

of  (1.16a,b)  for  ffi  that  the  sharp  coefficients  f  are  a  more  natural 

vehicle  for  our  algorithms  than  the  daggar  coefficients  f^.  Nevertheless, 

the  latter  are  of  algorithmic  and  theoretical  importance.  Indeed,  they 

00 

are  needed  for  a  final  inversion  f(i)  =  £  f^h  (t)  returning  to  the  func- 

.00 

t ion  sought . 

As  in  the  introductory  paper  [A],  the  general  setting  for  the  trans- 

00 

form  is  the  L_ (-“,“)  class  for  which  T  f+^  <  «,  so  that  (f  )°°  e  t~.  As 

2  L  n  ’  ^  n  -«>  2 

.00 

for  the  one-sided  functions,  however,  the  methodology  is  of  value  algo¬ 
rithmically  only  for  functions  sufficiently  smooth,  e.g.,  "rapidly  decreas 
ing”  (cf.  §2).  Summability  of  (f|)  is  of  special  importance.  We  note, 
therefore,  that 

(1.17)  (fV  £  t.  =>  (fV  e  l  =>  1  f#  =  0 

n  -oo  1  '  n  -oo  1  <-  n 


*S  ^  *f*  *1" 

as  seen  immediately  from  (1.14),  and  l  fn  =  fg  -  f^ 


(1-18)  (<)1  € 


i,  =>  r  =  l  f#  =  -  l  r  , 

1  n  L  m  L.  m 


n+1 


Moreover, 


#  t 

and  (1.18)  then  permits  one  to  go  from  (fn)  to  (fn) , 
and  (1.14)  that  for  all  (f^)  e  L ^ 


We  also  note  from  (1. 


(1.19) 


<=  < 


n+ 


f0+  +  f0 


(-n)-  ' 


n  >  0 
# 


,  n  =  0 
n  <  0  . 


4d) 
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Consequently  one  has 


Tr(u) 


-  I 


J  n 
f  u 
n 


0 

I 


(-n)- 


J 

f  u 
n+ 


n 


i.e. , 


(1.20)  T*(u)  =  T*+(u)  +  tJju'1)  . 


-9- 


§ 2 .  Rate  of  convergence  of  the  Laguerre  coefficients;  regularity  structure 
in  the  complex  s  and  u  planes 

A.  Structure  for  the  one-sided  case 

The  bilateral  Laguerre  transform  will  be  a  useful  tool  for  the  mechaniza- 

+  GO  ft  00 

tion  of  convolution  if  the  sequences  (fn)  ^  and  (f  )  ^  fall  off  rapidly  as 
in  [A].  The  definitions  employed  there  will  be  repeated  for  convenience. 


DEF  2.0a.  (f  )™  £  C*(N+)  n ^ |  f ^  |  •+  0,  n  -*■  00  for  all  non-negative  integers 

K. 

DEF  2.0b.  f(x)  £  C  +  (R+)  |T^(^-)^f  (t)  |  <  for  all  non-negative 

integers  q,p. 

00 

One  then  says  (f  is  a  rapidly  decreasing  sequence  and  f(r)  is  a 
rapidly  decreasing  function  on  R+. 

If,  in  particular,  f  (t)  is  "rapidly  decreasing",  then  as  in  [A],  Theorem 

6.6,  (f*)  and  (f*)g  will  be  rapidly  decreasing,  e.g.,  one  will  have  n^|f*|  -*•  0, 

n  °°  all  positive  integers.  K.  Similarly,  if  f  (-t)  is  rapidly  decreasing 
+  00  ft  00 

on  (0,®),  then  (f  ^  ^ ^  and  (f  n  j)0  will  also  be  rapidly  decreasing. 

A  systematic  development  of  the  extension  of  these  ideas  to  the  two- 
sided  setting  will  be  given  soon. 


P2.1 


+  00  00  \  ft  CO  CO  R  00 

(fj)0  £  C+(N+)  <£>  (fn)0  c  C+(N+)  <%  f(T)  £  C+ (R+) 
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Proof:  is  immediate  from  f*  =  f^  -  f^  ,  in  (1.14). 

-  n  n  n-1 


A-  nK|f+|=nK|-  l  l  m*+1|f"l  <  -±-  0  as  n  -  - 

'  n  nil  m  n  nil  m  n 


J,  _  1  r  K+li,#i  .  K 


for  all  K  . 


Restatement  of  Theorem  6.6  in  [A].  □ 


We  need  to  relate  our  smoothness  conditions  to  a  radius  of  convergence 
for  T^(u) .  We  define  the  following  set. 


DEF .  2.2 

n=0 

We  then  have 

P2 . 3a .  If  a  sequence  has  a  radius  of  convergence  larger  than  unity,  then 
the  sequence  is  rapidly  decreasing,  i.e., 


(2.1)  Gr  c  C+(N+)  ,  any  R  >  1 


Proof 


P2.3b. 


*K|fJl  =  (4)ifn!Rn  <  C  TT*  0  as  n  -  “  •  D 

R  R 

(<)o  £  GR  <=>  (0o  £  GR  >  for  any  R  >  1  . 


Proof :  The  direction  =>  follows  at  once  from  (1.14).  The  direction  <=  is 

an  immediate  consequence  of  T^.(u)  =  (1-u)  T^(u)  relating  the  generating 

ft  00  ft 

functions.  Since  (fn)Q  €  GR,  R  >  1  =>  T^.(u)  is  regular  for  0  <  |u|  <  R, 

the  point  u  =  1  lies  in  the  domain  of  regularity  of  T^(u).  Moreover, 

#  t 

TjO)  =  0.  It  then  follows  that  u  =  1  is  a  regular  point  for  T^(u),  etc.  □ 
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Remark  2.5c.  We  note  that  c”(N+)  c  ,  and  is  a  proper  subset  of  G^,  as 

Co  J 

may  be  seen  from  (fn)g  with  f  =  - j  •  We  also  note  that  the  union 

1  +  n 

of  the  sets  GR  with  R  >  1  is  a  proper  subset  of  C+(N+),  as  in  the  example 


$  ^  Tl 

We  will  see  that  over-convergence  of  T,(u)  =  £  f  u  ,  i.e.,  the  avail- 

t  0  n 

ability  of  a  radius  of  convergence  greater  than  one  is  associated  with  over 
convergence  of  the  Laplace  transforms  4>(s)  in  the  presence  of  a  simple  con¬ 
dition  at  infinity. 

P2.4.  Let  f (t)  e  L2( 0,°°)  with  <t>(s)»  Tf(u)  defined  for  Re(s)  >0,  |u|  <1. 

Then  for  the  analytic  continuation  one  has 

oo 

{<Ks)  =  /  e"STf(t)dT  is  regular  at  s  =  “  and  <j>0)  =  0)  <==> 

0 

#  # 

(T^(u)  is  regular  at  u  =  1  and  Tf(l)  =  0} 

Proof:  From  (1.12),  T*(u)  =  <K  j  j^)  and  <Hs)  =  Tf  ■  By  setting 

1  1  #  2-w  1 

w  =  <f> (— )  =  •  'ihus  4>(~)  is  regular  at  w  =  0  and  vanishes  there 

#  # 
if  and  only  if  T^(u)  i's  regular  at  u  =  1  and  T^Q)  =0.  □ 

The  situation  is  significantly  modified  when  4>(s)  is  not  regular  at 

infinity.  Some  simple  bilinear  mappings  from  the  s-plane  to  the  u-plane 

are  needed  to  understand  the  modified  behavior.  We  state  the  basic  results 

which  follow  from  standard  results  of  conformal  mapping  theory  [9]. 

P2.5. 

(A)  Let  s(u)  =  —  and  A  =  {u:  Rj  s  |u|  <  R^),  where  0  <  R^  <  1  £  Rj.  Then 
A  maps  into  s(A)  in  the  complex  s-plane,  the  complement  of  the  two  dis¬ 
joint  circles  possibly  tangent,  as  shown  in  Fig.  2.1. 
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(B)  Let  u(s)  =  ^  ~  j  and  B  =  {s:  Rj  s  |sj  <  R2},  where  0  i  Rj  ^  J  ~  ^2" 

Then  B  maps  into  u(B)  in  the  complex  u-plane,  as  shown  in  Fig.  2.2. 

(C)  Let  u(s)  -  -jg  and  D  =  {s:  6^  s  Re(s)  £  6^},  where  6^  £  0  £  6^.  Then 

u(D)  for  two  key  cases  is  as  shown  in  Fig.  2.3a,b. 

Two  key  theorems  for  the  one-sided  functions  may  now  be  stated. 

Theorem  2.6 

Let  T^(u)  =  <f>  (j  y^)  .  Then  the  following  (A)  and  (B)  are  equivalent. 

(A)  (f  )”  e  G_  for  some  R  >  1. 
v  '  v  n'O  R 

(B)  4>(s)  is  regular  in  {s:  6  <  Re(s)}  for  some  6  <  0.  <}>(s)  is  also  regular 

at  s  =  °°  and  4>(“)  =  0. 

Proof 

The  theorem  follows  from  P2.5  and  simple  argument  in  the  complex  plane. 
Details  are  omitted. □ 

If  the  regularity  condition  of  at  s  =  00  is  dropped  in  (B)  of  Theorem 

2.6,  T^(u)  is  no  longer  guaranteed  to  be  regular  at  u  =  1.  The  resulting 
regularity  structure  is  described  next.  The  proof  is  as  before  with  details 
omitted. 

Theorem  2 . 7 

Let  T*(u)  *  Then 

(1)  The  following  (A)  and  (B)  are  equivalent. 

#  * 

(A)  T^(u)  is  regular  inside  a  circle  c  containing  every  point  of  the 

1  1  * 

set  {u:  I u }  =  1,  u  *  1}  in  its  interior,  with  the  circle  c  tangent 
to  the  unit  circle  at  u  =  1 . 

(B)  $(s)  is  regular  in  =  {6:  6  <  Re(s)}  for  some  6  <  0. 
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(2)  Let  s  £  E.  be  such  that  s  =  S  ♦  iT.  Then  (B)  implies  lim  $(S  ♦  iT)  = 
#  T-±- 

lim  Tf(u)  =  0,  the  limit  being  taken  along  the  unit  circle, 
u-l 

Corollary  2.8 

If  4> Cs)  is  entire,  then  T^(u)  has  only  one  singularity  at  u  =  1. 

This  corollary  also  can  be  seen  directly  from  T^(u)  =  4> (y  y^-) . 

B.  Structure  for  the  two-sided  case 

Theorems  2.6  and  2.7  can  be  easily  extended  to  the  two-sided  case  in  a 
natural  way.  First  we  give  the  following  definition. 

DEF 

(2.2a)  f(T)  £  C*(Rj)  <=£  f+(T)  £  c"(R+),  f_(-x)  £  c“(R+) 

(2.2b)  (fn)l  e  C* (N)  iSt  (fn)J  €  c“(N+).  (f.n)“  £  c”(N+)  . 


Remark  2.9:  The  definition  implicit  in  (2.2a)  requires  no  continuity  at  t  =  0 

CO 

for  f(r).  The  asterisk  distinguishes  this  class  from  the  ordinary  class  C+(Rj) 
of  rapidly  decreasing  functions  which  are  differentiable  to  all  order  at 
every  t . 

These  definitions  with  P2.1  yield  the  following  immediately. 

P2.10 

(2.3)  (f*)".  ^  C* (N)  <=>  (/)"„  e  C*(N)  <=>  f(r)  £  C*(R)  . 


We  also  have  as  a  simple  extension  of  the  one-sided  case: 
P2.ll 


(fV  <  t,  =>  (f')“  € 
n  -«  1  v  n  -<*> 


l1  and  l 


f#  =  0 
n 
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Proof 


(f*)”^  e  Z ^  implies  that  (f*+)“  e  Z^  and  (f^_)Q  e  Z 1  and  that 


(f  )"  e  Z.  . 
v  n  -°°  1 


Since  T*(u)  =  T*+(u)  +  T*_(u-1)  from  (1.20), 


(2.4)  T  f#  =  l  f#  +  I  f*  =  0  .  □ 

L  n  Ln  n+  Ln  n- 
n=-oo  n=0  n=0 


The  extension  of  Theorems  2.6  and  2.7  to  the  two-sided  case  now  follows. 
Theorem  2.12 

Let  T*(u)  =  •  ^en  the  following  (A)  and  (B)  are  equivalent. 

(A)  T^-(u)  is  regular  in  AR  =  (u:  R  <  |u|  <  R+)  where  0  <  R _<  1  <  R+. 

(B)  <}>g (s)  is  regular  in  E  =  (s:  -A  <  Re(s)  <  A+)  where  A_  >  0,  A+  >  0,  as 
well  as  at  s  =  °°,  where  <j>g(“)  =  0. 

Proof 

Details  omitted. □ 

We  note  that  from  P2.3a,b  and  P2.10  that  either  (A)  or  (B)  of  Theorem 
2.12  implies  f(t)  €  C+(R).  The  two-sided  counterpart  of  Theorem  2.7  is  given 
next. 

Theorem  2,13 

Let  T*(u)  =  *B(~  ~).  Then 
(1)  the  following  (A)  and  (B)  are  equivalent. 

II  *  *  *  * 

(A)  T^(u)  is  regular  inside  the  region  C  *  ^2  ~  ^1  where: 

* 

(a)  C2  is  a  circle  containing  every  point  of  the  set 
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{u:  } u J  =  1,  u  *  1}  in  its  interior; 

*  * 

(b)  is  a  circle  where  -  {1}  is  contained  in  the 

interior  of  the  set  (u:  |u|  s  1}; 

★  ★ 

(c)  both  Cj  and  C  are  tangent  to  the  unit  circle  at  u  *  1. 

(B)  4> B Cs)  is  regular  in  E  =  {s:  -A  <  Re(s)  <  A+}  where  A  >  0  and 

A+  >  0. 

(2)  Let  s  =  S  +  iT  for  s  e  E.  Then  (B)  implies  that  lim  <J>_(S  +  iT)  = 

#  T-±- 

=  lim  Tf(u)  =  0,  the  limit  being  taken  along  the  unit  circle. 
u+1 

Proof : 

Omitted. □ 

C.  Rapidly  decreasing  functions  on  the  full  continuum 

P2.10  says  that  f(i)  c  C+(R)  if  and  only  if  (nKf^)  e  -fj  for  all  non- 

K  # 

negative  integer  K,  or  equivalently  (n  fn)  e  Z^  for  all  nonegative  integer 
K.  But  this  permits  the  discontinuity  of  (^-)rf(x)  at  t  =  0.  When  (j£-)rf(t) 
is  continuous  at  t  =  0  for  all  r  =  0,1,2,...,  one  finds  surprisingly  that 

X  ft 

all  the  moments  of  f1  and  f  vanish.  We  remind  the  reader  that  Z  (t)  satis- 

n  n  n 

fies  the  operator  relation  (cf.  Eq.  (6.2)  of  [A]) 

(2.5a)  L[£ri(t)]  =  (n  +  i)£n(T) 

where 

(2.5b)  L  -  [I  T  -  ^  . 

From  (1.4a)  and  (2.5)  we  have  L [h^ (t ) ] 

L[hn(t)]  -  -L[£_n.1(-T)U(-t)]  =  -  (-n 
Consequently  one  has 


=  (n  +  j)hn(t)  for  n  2  0.  For  n  <  0, 
'  1  +  iH.n.il-OlK-O  ■  (n  ♦  j)hn(0. 
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(2.6)  lr[hn(T))  =  (n  +  y)rhfl(T)  for  all  n  and  r  =  0,1,2,... 

We  also  verify  trivially  from  equation  (1.4a,b)  that 

(2.7)  0  =  h  (0+)  -  h  (0-)  =  1  for  all  n  . 

v  n  n  nv 

It  follows  that  for  f(i)  e  C.(R,), 

T  1 

00 

(2.8)  Lrf (t)  °q*  =  l  (n  +  |)V  ,  r  =  0,1,2,...  , 

n=-® 

0+  oo 

where  a(t)  Q  =  a(0+)  -  a(O-).  Legitimacy  is  assured  by  (n  f  )  m  e  for 

all  k  £  0  (cf.  Theorem  6.6  of  [A]),  the  uniform  boundedness  of  £^(t)  on 

Rj  (cf.  Lemma  6.2  of  [A])  and  the  dominated  convergence  theorem.  Similarly, 
•*•00  * 

when  (f  )  e  C  (N) ,  (2.8)  will  be  valid.  The  following  theorem  can  now 

n  »oo  t 

be  proven. 

Theorem  2.14 
Let 

00 

(2.9)  A  -  (f(T):  (f*)"„  e  C*(N)  ,  l  n  V  =  0  ,  k=  0,1,2,...}  . 

—  00 

Then 


C+(RX)  =  A  . 


We  note  that  A  is  a  proper  subset  of  C+(Rj)  which  is  equal  to  (f(x):  (f^)  ^  t 
* 

C+(N)}.  It  differs  in  that  its  elements  have  all  moments  equal  to  0.  The 

OD 

class  C  (R.)  contains  real  analytic  functions  of  great  interest  in  mathematical 

+  1  2 

—  x  /  2  —  i 

statistics  and  probability  theory,  e.g.,  e  and  (cosh  x)  .  The  vanish- 

•  k  -j.  “  k  ft 

ing  of  J  n  f  and  as  we  will  see  of  T  n  f  for  such  functions  is  of  corre- 
°  L'  n  **  n 

—  00  -oo 

sponding  interest.  Two  lemmas  will  be  employed  to  prove  the  theorem. 
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Lemma  2.15 


l  "k<-  ».  k-  <->  l  "k<”  »■  k  =  0,1,2,... 


Proof 


#  t 

(=>)  Immediate  from  f  =  f 

n  n 


(<=)  One  has 

(2.10) 


Tj(u) 


Tf(u) 


(1-u)  (1-u) 


k+1 


-  f^  ,  and  the  binomial  theorem, 
n- 1 


for  u  =  e10  ,  6  *  0 


The  limit  of  the  expression  on  the  right  exists  and  equals  0  as  u  -*■  1  by 

repeated  application  of  L'Hospital's  rule.  A  similar  application  of 

L'Hospital's  rule  to  the  expression  on  the  left,  then  gives  via  induction 
00 

on  K,  l  nKf*  =  0  for  all  K  as  required. □ 

-QO 

Lemma  2.16 


f(x)  e  A  => 


(a)  f(i)  exists  at  t  =  0 
<  .(b)  ^  f(t)  e  A 

t  (c)  (|T)Xf (x)  e  A 


Proof 

We  first  show  that  f(T)  e  A  implies  that  —■  f(t)  exists  at  t  =  0. 
This  is  easy  since  from  (2.5b)  and  (2.8) 


L[ f(x)] 


0+ 

0- 


=  I  (n 


I>f 


t 

n 


0  . 


Since  f(x)  e  C*^)  implies  that  f(x)  is  differentiable  in  (-®,  0)  and  (0,  «), 
differentiability  on  Rj  follows. 

We  must  now  show  (b)  of  the  Lemma,  i.e.,  that  f(i)  e  A.  One  has 
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for  continuously  differentiable  f(x)  that 

(2.1D  4f(T)];-K  +  0  * 


from  Theorem  4.4  of  Section  4.  From  (2.11)  and  the  binomial  theorem,  we  have 
(2.12)  I,.  nK[^  f(r)]J  -  i  j[  nV  ♦  I  ^  {(n-1)  ♦  l}1^  =  0  . 

Finally,  from  (2.12)  and  Lemma  2.15,  one  has 

£50 

l  ^[37  f  ("0  ]n  =  0  ,  K  =  0,1,2,...  proving  (b)  . 


Since  (c)  follows  immediately  by  induction,  the  Lemma  is  proven. □ 
Proof  of  Theorem  2.14 

10+ 


(=>)  f (r)  e  c”(Ri^  =>  LTf(r)  e  C^(R^)  so  that  lTf(t) 

art  co 

K-+ 


0- 


=  0 .  Thus 


l  (n  +  y)rf^  =  0  from  (2.8),  whence  \  n  f^  =  0  for  all  K. 


(<=)  It  is  easy  to  see  that 


f(t)  e  C+(Rj) 


(2.13)  f(x)  e  C+(Rj)  <=> 


0+ 


=  0,  r  =  0,1,2,... 


Consequently  we  need  only  verify  that  every  derivative  of  f(t)  Is  continuous 
at  x  =  0.  The  theorem  therefore  follows  from  Lemma  2.16  and  (2.13).D 
Remark  2.17:  It  is  natural  to  ask  how,  in  the  setting  of  Theorem  2.13,  the 
continuity  of  f(x)  and  all  its  derivatives  at  x  =  0  modifies  the  behavior  of 
$g(s)  in  its  convergence  strip.  The  answer  is  given  by  Theorem  (a)  of  §2.2 
in  Dym  and  McKean  [1]  which  says  that  the  Fourier  transform  of  a  function 

ao  cd  ♦ 

in  C+(Rj)  is  also  in  C+(Rp.  For  us  this  means  that  $  (u)  =  ^g(iu)  is  in 
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c”(Ri) ,  i.e.,  that  (u)  0  as  |u|  ®.  It  is  clear  from 

ao 

j  e~STe~aTf Ct) dx  =  <f>g(s+a)  that  every  function  +  iu)  for  a  in  the 

-00 

interval  of  absolute  convergence  of  the  bilateral  Laplace  transform  will 
also  be  rapidly  decreasing. 

Remark  2.18:  A  discussion  of  rate  of  convergence  of  the  Laguerre  transform 
coefficients  for  functions  having  only  finite  number  of  drivatives  may  be 
developed  along  the  line  of  Section  6  of  [A].  This  development  will  be 


omitted. 
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3 .  The  extent  of  the  Laguerre  transform  coefficients 

The  accuracy  of  the  algorithm  based  on  the  Laguerre  transform  depends 

on  one's  ability  to  represent  the  functions  present  compactly,  i.e.,  with 
,  K 

■J*  4- 

vectors  (f  )  v  of  reasonable  length.  Correspondingly,  some  measure  of 
n  -  K+  ' 

extent  of  the  sequence  (fT)  may  be  useful  and  feeling  for  the  relation- 
ships  between  the  extent  of  the  approximating  sequence  and  that  of  the 
function  helpful.  In  general,  there  are  inverse  relationships  between 
these  two  extents  as  in  the  Heisenberg  inequality  of  Fourier  transforms 
(cf.  Dym  and  McKean  fl],  §2.8).  In  what  follows,  we  discuss  only  rapidly 
decreasing  functions  f(t)  in  either  C.(R)  or  C°°(R)  .  This  condition  may 

*♦  ■f 

be  weakened,  as  mentioned  in  Remark  2.18. 

4-  00 

We  first  exhibit  such  inequality  for  the  dagger  transform  (f  )  .To 

n  ~ 00 

do  so,  we  require  the  following  notation.  Let  be  the  extent  of  the 
t  00 

dagger  transform  (f^)  ^  defined  by 

(3.1)  N l*  l  )n  +  1/2 1  f f/[  ff  . 

-00  -oo 

Similarly,  let  Tj  be  the  extent  of  the  original  function  f(-r),  i.e., 

00  00 

(3.2)  T2  =  /  |t|f2(T)dT/f  f2(i)dt  . 

-00  _00 

Finally,  let 

00  00 

(3.3)  Bj  =  /  )t Jf’ (T)2dr//  f2(i)di  . 

-00  -oo 

Then  one  has  the  following  theorem. 
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* 


s 


Theorem  3 . 1 


Let  f (t)  e  C  (R) .  Then 

T 


»)  -  j  f,  *  ; 


(2)  Rj  ,  1(1,  .i-)*! 
*  1 


Proof 


From  (2.6),  L [f (t ) ]  =  J  (n  +  j)hn(j)  so  that 

-00 

00  00  .2 
/  L[f(T)]f(r)dT  =  l  (n  +  i)fT  >  0 
0+  n=0  2  n 

and 


0-  -1  2 

/  L[f(r)]f(T)dT  =  l  (n  +  y)ff  <  0 
-°°  n=-°° 


We  note  from  integration  by  parts  that 

00  CO 

/  L[f(T)]f(T)dx  =  /  {j  Tf(T)  -  4-  Tf' (x)}f(T)dT 

0+  •  0+  dT 


1  -  O  w 

~  J  !  Ct)dT  +  /  Tf  (t) 2dl 
0+  0+ 


Similarly, 


0-  0- 

-/  L[f(T)]f(T)dT  =  /  {j(-T)f(T)  +  g-  Tf'(T)}f(t)dT 


1  0-  0-  , 

*4  /  Mf2(T)dT  +  /  I  T  I  f '  (T)  2dx 

-oo 


AO-A085  100  CENTER  FOR  NAVAL  ANALYSES  ALEXANDRIA  VA  OPERATIONS  EV— ETC 
THE  LA6UERRE  TRANSFORM. (U) 

MAY  80  J  KEILSON*  W  NUNN»  U  SUMITA 
UNCLASSIFIED  0E6-PP-284 


F/6  12/1 

ML 
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Hence , 


i  i"  *  ii<2  ■  !  (» *  X2  -  i1  <■  *  X2 

-oo  n=0  -® 


00  CD 

=  J  /  |  T  |  f 2  (t)  dt  +  /  I  T  I  f  (t)  2dt 


00  2  t 

Dividing  both  sides  by  /  f  (t)dt  =  J  f  ,  one  has 


(3.4)  nJ  =  I  T*  +  Bj  , 


proving  part  (1).  To  prove  part  (2),  we  first  note  the  identity 


00  °0  00 

(3.5)  /  Tf'(T)f(T)dt  «  I  /  T  f2(T)dT  =  -  i  /  f2(T)dT 

-00 


By  Schwarz's  inequality, 


jJ  f2(T)dt  <  /  | t | |f' (t) { |f(x) |dx 


so  that 

(3.6)  yStj®! 


Substituting  (3.6)  into  (3.4),  one  has  the  desired  result .0 
Remark  3 . 2 

We  see  from  Theorem  3.1  that  when  Tj  becomes  infinite  or  when  goes 
to  zero,  Nj  becomes  infinite.  The  methodology  therefore  cannot  tolerate 
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functions  f  (x)  too  closely  concentrated  at  zero  or  functions  f(x)  too 
great  in  extent. 

The  inequality  of  Theorem  3.1  has  a  counterpart  for  the  sharp  trans- 

^  QO 

form  (fn)  ^  of  equal  interest.  As  we  will  see,  the  sharp  transform  is 
less  sensitive  to  concentration  of  f(x)  at  t  =  0,  a  useful  and  important 
advantage.  The  following  identity  underlies  the  inequality  for  the  sharp 
transform. 

Theorem  3.3 

Let  f(x)  e  C+(R).  Then 

l  lnlfn  =  /  Mf2(r)dT  . 

.00  .00 

Proof 

Let  v(x)  =  hg (t) *f  (x)  so  that  f(x)  =  v  (x)  +  j  v(x).  Then 
00  00 

(3.7)  /  |x|f2(x)dx  =  /  |x|{v'(x)  +  \  v(x))2dx 

.00  -  CD 

00  *  00  oo  ^ 

=  /  | T I v  (x)2dx  *  J  f  |t|v(x)  d X  +  /  I x I  V  (x)v(x)dx  . 


Since  v(x)  is  continuous  at  x  =  0,  as  we  will  see  in  Theorem  4.2  of 
t  # 

Section  4,  one  has  v  =  f  for  all  n.  Hence  from  Theorem  3.1,  the  sum  of  the 

n  n  n  ~ 

*  i  +2  00  i  #2 

first  two  terms  in  (3.7)  is  equal  to  J  |n  +  y|v  =  l  |n  +  y|f  . 

.00  .00 

other  hand,  the  third  term  becomes,  by  an  analogy  to  (3.5), 


/  Mv  (x)v(x)dx  *  /  xv  (x)v(x)dx  -  j  xv  (x)v(x)dx 
.00  0  .00 


0 


-  j  I  V2(x)dx  +  \  I  v2 (x)dx 


On  the 
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•1  J- 


1  P  t  1  r  t 

=  '  2  l  Vn  +  2  l  \ 

U  -OO 

oo  2  oo  2  oo  2  -1  2 

Since  I  |n  +  -||f*  =  l  |n|f*  +  ^  I  ~  j  l  v*  ,  one  has  from  (3.7)  that 


/  Mf2(t)  =  l  |n|fn  . 


D 


Theorem  3.3  now  leads  to  the  counterpart  of  Theorem  3.1  for  the  sharp 

transform  (f#)“  .  As  in  (3.1)  for  the  dagger  transform,  we  define  for  the 
n  -® 

sharp  transform 


(3.8) 


Theorem  3.4 

Let  f (t)  e  C* (R) .  Then 


(1)  N*  =  pTj  where  p  =  l  f *  /f  f* 

-OO  -00 

(2)  N*  >  |  T j  . 


Proof 

Part  (1)  is  immediate  from  Theorem  3.3.  For  part  (2),  one  has  from 
Schwarz's  inequality 


-00  -oo 


*  4 1 


n 


so  that  p 


The  result  then  follows. D 


Remark  3.5 


We  see  from  Theorem  3.4  that  when  becomes  infinite,  so  does  Nj.  But 

in  contrast  to  dagger  transform,  concentration  of  a  function  f(T)  near  zero 

-# 

does  not  necessarily  imply  that  Nj  becomes  infinite.  This  point  will  be 
seen  vividly  through  the  following  example.  Let 

f(r)  =  {ee"6Tu(T)}  +  {ee0Tu(-T)}  ,  e  >  o  . 


We  note  that  the  larger  6  is,  the  more  f(t)  is  concentrated  near  zero  or, 
equivalently,  the  smaller  6  is,  the  greater  the  extent  of  f(t)  is.  One  finds 
readily  that 


*  +2 

(3.9)  l  f ;  *  0  ; 

_oo 

oo  2  oo 

(3.10)  l  jn|f*  =  /  | r | f 2 (t) dt  =  ~  ,  for  all  6  >  0  ; 

— oo  .oo 

00  2 

(3.11)  l  f*  *  40 (0+1)/ (26+1) 2  . 

.00 

From  (3.9),  (3.10)  and  (3.11),  one  has 


(3.12) 


(3.13) 


(3.14) 


N#  -  1  (29+l)2 
W1  "  2  40(0+1) 


P 


(26+1) 2 
4(6+1) 


We  see  easily  that  p  is  a  monotone  increasing  function  of  8.  A.t  6  *  0, 
p  *  1/4  and  therefore  p  2  1/4,  as  expected.  We  note  that  the  inequality  (2) 
of  Theorem  3.4  is  thus  sharp.  As  0  +  0,  both  and  Nj  go  to  +«®.  But  as 
0  ■+  +<*>,  fjj  -*■  1/2  while  Nj  -+  +»  since  -+  0. 
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+  # 

The  relation  v  =  f  with  v(x)  =  hn(x)*f(x)  played  a  key  role  in  the 

n  n  u  “  #2 

proof  of  Theorem  3.3.  This  also  enables  us  to  evaluate  £  f  in  a  closed 

.00 

f  #  1 

form.  The  relation  g  =  nf  with  g(x)  =  L[f(x)]  -  f  (r)  given  in  Theorem 

00  2 
2  a 

4.6  of  Section  4  leads  to  a  closed  form  of  )  nr 

L  n 
.00 

Theorem  3.6 


Let  f(x)  e  C,(R).  Then 

Y 


U) 


00  °°  .ill 

\  f n  *  /  e  'i'T'r^.(T)dT  where  r^.(x)  =  f(x)*f(-x) 


00 

Furthermore,  if  f(x)  is  in  C+(R), 

(2)  l  n2f*  =  j  S  T2f2(T)dt  +  /  r2f  (r)2dr  . 

.00  n  -»  -°° 

Proof 

Let  v(t)  =  h0(i)*f(T)  so  that  =  f*  as  before.  Then 

«  #2  «  2  •  , 
l  fn  =  l  'vn  =  /  {h0(T)*f(Tj}2dT  . 

.00  -CO  -00 

00  oo 

Let  4>(u)  =  /  e1UTf(r)dT  and  h(u)  =  /  e*UTh_(x)dT.  By  the  Parseval  identity, 

—  OO  _oo 

one  has 


l  =I?/  |h(u)|2|<^(u)|2du 

.00  —OO 


=  /  {h0(T)*h0(-x)}{f(x)*fC-T)}dT 

.oo 


e'^l 


rf(t)dx 
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proving  (1).  To  prove  part  (2),  we  let  g(x)  =  L[f(-r)  -  f  (t)] .  Then,  as 

■j*  |  qo 

given  in  Theorem  4.6,  one  has  g  .=  nf  for  f(x)  e  C,(R).  Therefore,  again 

n  n  ▼ 

by  the  Parseval  identity, 


00  2  00 

l  nV  =  /  (L[f(x)]  -  f’(x))2dx 


=  +  u2)|^- <f(u)|2du 


du 


*  j  f  t 2f 2 (x) dx  +  /  {f (t)  +  Tf  (t) }2dx 

—  OO  _  CD 

oo  oo 

=  J  /  T2r(T)dT  +  /  T2f'(T)2dT  .  □ 

-00  —  oo 


Theorem  3.6  gives  relations  for  the  second  moments  of  the  extent  of  the 
sharp  transform  and  the  original  function.  Let 


(3.15) 

(3.16) 

(3.17) 


2 


T 


2 


B 


2 


-QO  — OO 


00  00 

/  x2f(x)dx//  f2 (x)dx  ; 

—  00  -.CO 


/  X2f'(T)2dT //  f2(T)dT 
-00  -00 


One  can  easily  show  that 

2 

(3.18a)  fi*  *  Nj 

(3.18b)  T2  a  f2  . 


Theorem  3.7 


Let  f(-c)  e  C°(R). 


(1) 

N#  = 
-# 

"<{  V  S2) 

(2) 

N2 

(|t2  .  B2)/f1 

(3) 

Tti 

s  fjj  s  (i  T2  .  82)/tj 

(4) 

S' 

( j  r  ♦  i)N*  where  r  =  T^/Tj 

Proof 

-# 

Part  (1)  is  immediate  from  Theorem  3.6.  From  Theorem  3.4,  =  pT^  and 

part  (2)  follows.  From  (3.18a), 

Ni  *dr-  ^t2  +  B2^f  • 

Ni 

1  -  -# 

Tj  -  is  proven  in  Theorem  3.4,  and  therefore  part  (3)  is  shown.  To  prove 
part  (4),  we  first  show  the  inequality 

(3.19)  T2  <  T2B2  . 

Let  sign  (t)  =  1  for  t  >  0  and  sign  (t)  =  -1  for  x  <  0.  Then 

■»  CO  0 

/  Sign(T)f ' (x)f(x)T2dT  =  /  X2f ' (x)f(x)dx  -  /  X2f' (T)f(T)dT  . 

0  -00 

Integrating  by  parts,  one  has 


/  T2f'(T)f(T)dT  =  J  /  T2  ~  f2(x)dx  =  -  /  xf2(x)dx  . 

0  o  aT  0 
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0  -  ,  0  0  7 
Similarly,  -  /  t  f  (T)f(r)dT  =  /  Tf^(i)dx  =  -  J  |i|r  (t)dT,  so  that 


00  CO 

,2 


j  ] t | f  Ct) dx  =  /  sign(T)f  (t) f ( T) t  dx  . 


Hence,  by  Schwarz's  inequality, 


J  |t|fZ(T)dT  <  /  /  T2f'(T)2dT  //  T^f(x)^dx  . 


2„,_,2 


2 

Dividing  both  sides  by  /  f  (x)dx  and  then  squaring,  one  has  (3.19).  Part  (4) 

.00 

follows  from  Part  (2)  and  (3.19).D 
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4.  Operational  properties  of  the  bilateral  Laguerre  transform 

In  principle,  operational  rules  for  the  bilateral  Laguerre  transform 

can  be  obtained  from  the  rules  for  one-sided  support  (cf.  [A],  Appendix  C) . 

One  works  on  f  (t)  =  f(-r)U(T)  and  f  (t)  =  f(T)U(-r)  separately  and  then 

00 

patches  the  results  together  carefully.  In  particular,  /  f(y)dy  may  be 

00  CO  f 

evaluated  in  that  way  only,  when  /  f(T)dx  %  0  and  /  f(y)dy  $  ") . 

-®  T 

The  operational  rules  for  the  bilateral  Laguerre  transform  are  of  separate 

interest,  however,  as  we  will  see.  They  contain  those  for  one-sided  support 

as  a  special  case.  We  also  derive  some  new  results  which  are  of  theoretical 

and  practical  value.  As  in  Section  3,  we  discuss  only  rapidly  decreasing 
★ 

functions  f(x)  in  C+(R).  All  the  results  are  summarized  concisely  in  a  table 
in  Appendix  A. 

We  first  extend  (1.15)  and  (1.16). 

P4.1 

★ 

Let  f(t)  e  C+.  Then 


g(i)  =  f(-T>  <  =  >  g*  =  f# 


gn  '  -f- (n+1) 


Proof 


00 

Since  Yfi(s)  =  tg(-s)  -  -  l 

—  CO 


f’(^i)n 
n  s-V 


one  has 


■r#c.i  rl  1+U-.  n 

V  2  1 -u  ^ 


no 

h  I 


,.+  -n 
f  u 
n 


T^u"1) 


Similarly,  T*(u)  =  ~  T*(u)  =  -  i^fu'1).  □ 

When  a  function  is  convolved  with  one  of  the  building  block  functions 

hm(T),  the  dagger  coefficients  of  the  resulting  function  become  the  sharp 

coefficients  of  the  original  function  shifted  by  -m.  We  see  this  next.  For 
00 

a  sequence  (a  )  ,  the  first  difference  will  be  denoted  by  Aa  =  a  -  a  ,, 

^  v  n  n  n  n-l 


I 

.« 
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and  the  second  difference  by  A  a  =Aa  -  Aa  =  a  -  2a  +  a 

1  n  n  n-1  n  n-1  n-2 

Theorem  4.2 

* 

Let  f(t)  e  C^(R).  For  all  integers  m,  one  has 

g (t)  =  h  (x ) *f (t)  <=>  g  -  Af#  ;  g+  =  f# 

6  J  nr  5n  n-m  ’  sn  n-m 


3n 


Proof 


It  is  easy  to  see  that  (u)  =  (l-u)um  for  all  integers  m.  Then  T#(y) 

m  & 


(l-u)umT^.(u)  and  the  result  follows .□ 


Of  particular  interest  is  the  case  m  =  0,  which  provides  a  closed  form 

00  2 
r  # 

of  fn  »  as  given  in  Theorem  3.6.  The  next  corollary  is  immediate  from 

—  CO 

Theorem  4.2. 

Corollary  4.3 

Let  f (x)  e  C* (R) .  Ther 

y 

(1)  8«  -  IM  -  h0(,)*fM  <■>  8*  =  /  j  ;  gj  . 

(2)  g(,)  -  f(T)  .  h.jWfW  <=>  £  .  ;  g l  -  ffntl  . 

The  bilateral  Laguerre  transform  of  a  derivative  of  a  function  is  a 
straightforward  extension  of  the  one-sided  case  with  slight  modification. 
Theorem  4.4 


Let  f(x)  e  C+(R).  Then 


g(T)  =  ~  f(T)  <=> 


.*1  -X  *  -  s„.o{£(0*>  - 


f7<-  l  <■  »*« 

m=n 


'4  =  s 


1  f+  +  y  f+ 

2  n  m  ’ 


n  <  0 
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Proof 

This  is  immediate  from  the  corresponding  rule  for  one-sided  support 
([A],  Appendix  C) .□ 

+  # 

We  have  seen  in  Theorem  4.2  that  g^  =  f^  with  g(x)  =  hp(T)*f(t).  The 

next  theorem  gives  the  reversed  relation  of  this,  i.e.,  functions  g(x)  and 
p  + 

f(x)  related  by  g  =  f  . 
v  r.  n 

Theorem  4 . 5 

Let  f (x )  c  C^(R).  Then 


g (T)  =  f(x)  <=> 


A  '  fl  -  Voff(0*>  -  f'°-» 
x  ■  -  jjK  -  Vo(f(o*>  -  ■ 


Proof 

From  Theorem  4.4, 


=  <-i>  -  s„,o(«0*>  -  f<°-» 

■  <  -  Vo'*'0*5  -  •  D 

We  recall  the  operator  L  -  x  -  x  ^-]  given  in  (2.5b).  For  functions 

★ 

f(x)  e  C,(R),  one  has  immediately  from  (2.6)  that 

r 


(4.1)  g(T)  =  L(r)[f(x)]  <=>  g*  =  M(n  ♦  l/2)rf*}  ;  g^  =  (n  ♦  l/2)rf* 


for  r  =  0,1,2,...  .  This  leads  to  the  following  theorem  which  is  of  theoretical 

00  2 

,  ...  .  ,  .  r  2  J  given  an  Theorem  3.6. 

value,  providing  a  closed  form  of  \  n  f  6 

Theorem  4,6  ~°° 

Let  f(x)  €  C* (R) .  Then 
▼ 


g(T)  .  L(f(l)]  -  ^  f(T)  <»» 


4  -  "<  *  -  ft°-» 


gf--  E,l-f.  *  Vo,f(0*>  - 

6n  n+1  ’ 


Proof 

From  (4.1)  with  r  =  1  and  Theorem  4.4, 

-  <”  -  1/2K-1  -  K  *  fn-l>  *  ’  f‘°-» 

=  nAf^  ♦  «n(0{f(0+)  -  f(D-)>  =  nf*  ♦  «n>0{f(0+)  ’  f(°-^  •  0 

The  bilateral  Laguerre  transform  of  tfCt)  is  exactly  the  same  as  for  the 
one-sided  case.  The  proof  is  straightforward  and  omitted. 

Theorem  4 . 7 

Let  f (t )  e  C* (R) .  Then 

g(T)  =  Tf(r)  <=>  g*  =  -  A2[(n+1)/+1]  ;  =  -  A[(n+l)f*+1)  . 

From  Appendix  A  of  [A] ,  one  finds  that 


(4.2)  I  hn(^un  =  exp[- 
-00 

with  the  understanding  chat  Eq.  (4.2)  is  valid  for  |u|  <1  when  x  >  0  and 
for  |u|  >  1  when  t  <  0.  This  leads  to  the  operational  rule  for  shifting 
given  in  Theorem  4.9.  A  preliminary  remark  is  needed. 

Remark  4 . 8 

Even  though  (h  (T))  {  one  sees  easily  that 

n  — oo  ^ 


T  1+U! 


(4.3a)  (Ahn(T))_K  e  t2 
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00 

(4.3b)  l  (Ahn(T))2  =  1 

_oo 

(4.3c)  |Ahn(T)|  <  1  for  all  integers  n,  T  *■  0  . 


These  statements  follow  from 


(4.4) 


[h0(T  -  T)]n  =  Ahn(T)  , 


as  the  reader  will  verify,  and 


(4-5) 


/  h2(t  -  T)dx  =  /  h2(T)di 


Theorem  4 . 9 


Let  f(t)  e  C+(R).  Then 


Proof 


g(T)  =  f(T  -  T)  <=>  g  =  l  r  Ah  (T) 

n  u  n-m  nr  J 


Clearly  y  (s)  =  e'sT*_(s)  so  that 


00 

r  # 

=  *  l  gn 


n+1 


T^(u)  =  d-iO^expC-^jT^) 


(1-u)  l  hn(T)unT^(u) 


=  I  (Ahn(T))un.T^(u)  , 


and  the  result  follows.  □ 
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5.  Applications 


In  applied  probability  theory  and  statistics,  the  prevalence  of  a 
k 

random  variable  S,  =  £  £.,  where  £.  are  i.i.d.  random  variables,  is 

j  =  l  J  J 

familiar.  One  also  often  encounters  the  ergodic  green  density  g(t)  = 

00 

no 

6(x)  +  l  av  } (t)  in  expressions  for  the  ergodic  distributions  of  homo- 
geneous  processes  modified  by  boundaries  [4],  [5],  Here  a  (r)  denotes 
the  k-fold  convolution  of  a  two-sided  p.d.f.  a(x)  with  itself.  The 
bilateral  Laguerre  transform  method  provides  an  algorithmic  basis  for 
calculating  multiple  convolutions  of  functions  with  two-sided  support  and 
thus  enables  one  to  evaluate  the  distribution  of  or  the  ergodic  green 
density  g(x)  numerically.  This  in  turn  is  the  key  to  numerical  evaluation 
of  a  variety  of  results  in  applied  probability  that  have  been  available 
only  formally. 

K 

In  the  expression  S..  =  ]  £.,  when  £.  are  centered,  the  central 

j  =  l  •* 

limit  theorem  says  that  SK  with  suitable  normalization  converges  to  a 
standard  normal  random  variable  in  distribution.  In  some  applications, 
the  question  of  importance  is  "How  fast  is  the  convergence  to  normality?". 
This  can  be  also  answered  by  the  method  directly  via  computation  of  the 
required  multiple  convolutions. 


In  this  section  we  discuss  the  algorithmic  procedure  for  the  calcula¬ 
tion  of  multiple  two-sided  convolutions  and  present  a  numerical  example. 

This  procedure  is  then  used  to  quantify  the  Lindley  process.  For  the  pro¬ 
cess  in  queueing  contexts,  the  ergodic  waiting  time  distributions  for  M/G/l 
and  G/G/l  systems  are  evaluated  thereby.  All  computations  were  done  on  a 
DEC  10  computer  in  a  time-sharing  mode  using  APL  as  the  programming  language 
The  DEC  10  APL  implementation  is  a  double  precision  system,  which  uses  a 
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precision  of  18  decimal  digits.  Relevant  formulas  were  usually  coded  in 

a  straightforward  way,  with  no  attempt  made  to  optimize  the  subroutines 

for  speed  or  accuracy.  In  spite  of  this,  the  results  displayed  here  were 

typically  obtained  with  CPU  times  in  seconds  to  at  most  a  few  minutes  and 

with  no  evidence  o.f  numerical  problems. 

Example  1:  Multiple  two-sided  convolutions 

Two-sided  convolutions  are  usually  handled  by  Fourier  transform 

methods.  There  does  not  appear  to  be  any  significant  literature  devoted 

to  numerical  evaluation  of  two-sided  convolutions.  When  the  function  being 

* 

convolved  is  rapidly  decreasing,  i.e.,  in  the  class  C+(R)  defined  in  Section 

3,  our  method  yields  a  fast  and  accurate  algorithm.  To  generate  multiple 

#  #  # 

convolution,  we  make  use  of  (1.13)  which  states  that  T^^(u)  =  T^.(u)Tg(u). 
Hence  (f*)  =  (a^)  ^  where  f(t)  =  a^(t),  and  one  has  the  following 
algorithm: 

I.  Representing  a(t)  as  a(t)  =  a+(t)  +  a  (t)  as  in  Section  1,  generate 

.  t  N  f  N 

or  store  in  the  computer  the  coefficients  (an+)0  and  (a^  )Q,  usually 

obtained  analytically  as  in  [A], 

-j-  -j.  #  £ 

II.  Convert  (a^+)  and  (a^  )  to  (an+)  and  (an  ),  respectively.  This 
corresponds  to  a  simple  differencing  operation  on  each  set  of  coefficients, 

since  T#(u)  =  (l-u)Tf(u) . 

#  N  #  #  #  #  # 

III.  Obtain  (a^)  ^  by  setting  a^  =  a^+  for  n  >  0,  a^  =  aQ+  +  aQ  ,  and 

#  # 

an  =  a(  n)  ^°r  n  <  0>  as  (1-19)- 

IV.  Perform  m-fold  discrete  convolution  on  (an)  Retain  only  2N+1 

#  N 

terms,  centered  at  n  =  0,  in  each  convolution.  The  result  is  (f  )  ... 

n  -N 
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V.  Convert  (f*)  to  (f^)  by  (1.18).  This  is  the  inverse  of  the  differ¬ 
encing  operation,  i.e.,  summation. 

ft  N  + 

VI.  One  has  fT  =  fT  for  n  2  0.  Hence,  sum  the  series  I  ft  (t)  to 

n+  n  n=0 

get  f+(t).  Similarly,  one  has  f^  =  -t  n  ^  for  n  s  0,  by  virtue  of  (1.4d). 

Hence,  sum  the  series  -  £  ^-n-l^n^_t^  tQ  ®et 

n=0 

This  procedure  is  only  slightly  more  complicated  than  that  for  the  one¬ 
sided  convolutions  given  in  [A] ,  but  the  CPU  requirements  for  the  new  pro¬ 
cedure  is  increased  approximately  by  a  factor  of  2.  When  a  function  a(t)  is 
symmetric  about  zero  so  that  a+(t)  =  a  (-t),  however,  one  can  evaluate 
multiple  convolutions  by  working  only  on  one  side,  as  the  reader  will  verify. 
Then  the  CPU  requirements  will  be  almost  the  same  as  for  the  one-sided  con¬ 
volutions. 

We  illustrate  the  procedure  with  a(t)  =  (2e  2tU(t) }*{etU(-t) } .  (a*+) 

# 

and  (a^  )  are  available  analytically  from  [A],  Appendix  B.  Fig.  5.1  displays 
a^(t)  for  1  $  it  <  7  and  -10  <  t  s  10.  Two  hundred  terms  each  for  a+(t)  and 
a  (t)  provided  an  accuracy  of  at  least  12  digits  after  the  decimal  point 
uniformly  over  the  interval  -10  s  t  s  10. 

The  procedure  is  also  applied  in  [11]  to  numerical  evaluation  of  multiple 
convolutions  of  the  Logistic  variate  with  the  p.d.f.  a(t)  =  e  V(1  +  e  t)2, 

_oo  <  t  <  ®.  The  Laguerre  transform  approach  appeared  to  be  more  systematic 
and  efficient  than  an  existing  method  [2]  which  relied  on  the  special  analytic 
feature  of  the  Logistic  variate.  The  Laguerre  method  also  evaluated  multiple 
convolutions  of  the  folded  Logistic  variate  as  a  by-product,  which  the  exist¬ 
ing  method  could  not  provide. 


a(x)  =  {2e'2Tu(T)}*{eTU(-r)} 


k  =  1 


0  2  4  6  6  10 


■sided  convolutions  a1  1  (t) 
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Example  2:  The  Lindley  process 

Consider  the  Lindley  process  [6],  i.e.,  a  random  walk  on  the  nonnegative 
half  line  defined  by 

(5.1)  Wk+1  =  [Wk  +  5k+1)  ,  k  =  0,1,2,... 

where  [X]  =  max{0,X},  is  a  random  variable  with  known  distribution  and 

5k  are  i.i.d.  and  independent  of  Wq.  We  assume  that  the  common  distribution 
of  £k  is  absolutely  continuous  with  density  function  a(x)  e  °°) .  We 

see  that  the  distribution  of  Wk+^  is  the  convolution  of  the  distributions  of 
W^  and  Ck+1>  modified  by  the  [  ]  operation.  The  effect  of  this  operation 
is  to  "sweep  up"  the  probability  to  the  left  of  the  origin  into  a  mass  point 
at  the  origin.  Let 


(5.2a)  Ek  =  P[Wk  =  0] 


x  <  0 


(5.2b)  Fk(x)  =  P[WR  s  x]  =  < 


[Ek  +  /fk(y)dy  x  >  0 


If  we  define  =  Wk  ♦  £k+j,  then  the  density  of  W^+1  is  given  by 


k+1 


(5.3)  C  (x)  =  Eka(x)  +  fk(x)*a(x) 


a  u 

Hence  if  Ek  and  the  Laguerre  sharp  coefficients  (a*)"^  and  (f*(k))p  are 
known,  one  has 


(5.4) 


_  # 
E,  a 
k  n 


m=0 


a  f  (k) 
n-m  nr  ’ 


-«•  <  n  <  • 
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00 

Obviously,  f^+1(x)  -  f^+1(x)U(x).  Hence,  since  £  f*  =  0  (cf.  (1.17)), 

—  CO 

(5.S)  fj(k+l)  =  l  fJJ#(k+l)  ;  f*(k+l)  =  fJJ#(k+l)  ,  n  *  1  . 


Finally,  E.  .  =  1  -  /  f.  ?(y)dy  =  1  -  2  J  (-l)nr  (k+1)  since  /  l (x)dx 
*  1  0  K  1  n=0  n  0  n 

(-l)n2  (cf .  [A] ,  Appendix  A) .  This  leads  to 


(5.6)  E.  .  =  1  +  2  2  f-  .  (k+1)  . 


Hence,  if  EQ,  fp(x)  and  a(x)  are  known,  we  have  an  iterative  schema  for  calcu¬ 
lating  E^  and  f^(x)  via  (5.4),  (5.5)  and  (5.6). 

We  illustrate  the  procedure  in  a  queueing  context.  A  stream  of  customers 
arrive  at  a  single-server  queue  at  a  sequence  of  arrival  epochs  Tq,  t^,  etc., 
the  k-th  customer  arriving  at  epoch  The  interarrival  times  T^  = 

are  i.i.d.  with  common  distribution  T(x).  The  service  times  required  by 
the  k-th  customer  form  a  separate  sequence  of  i.i.d.  random  variables  with 
common  distribution  S(x),  where  and  T^  are  independent.  If  is  the 
time  the  k-th  customer  must  wait  in  queue  for  service,  then  one  has  [6] 

(s.7)  Vj  ‘  lwk 4  W*  '  H.i  ■  sk  -  Tk  •  k  - 

This  process  is  called  the  Lindley  waiting-time  process.  Clearly,  the  common 
density  of  is  given  by 

(5.8)  a(x)  -  (s(x)U(x) }*{t(-x)U(-x) } 


where  s(x)  and  t(x)  are  the  density  functions  for  the  service  time  and 


interarrival  time  variables,  respectively.  When  the  system  starts  with 
an  empty  queue,  one  has  Eq  =  1  and  f^(x)  =  0. 

The  iterative  schema  (5.4),  (5.5)  and  (5.6)  was  first  tested  in  the 
M/M/1  setting  for  which  the  limiting  distribution  of  as  k  goes  to 
infinity  is  analytically  available  [7].  For  the  traffic  intensity  p  =  0.5, 
an  accuracy  of  7  digits  after  the  decimal  point  was  attained  at  k  =  50, 
while  for  p  =  0.67,  k  =  128  was  needed  to  attain  the  same  accuracy  level. 

The  procedure  was  then  applied  to  the  M/G/l  queueing  system  with  s(x)  = 

_  x  /2  1 

xe  '  ,  the  Rayleigh  distribution  and  t(x)  =  e~  .  The  survival  functions 

of  for  various  values  of  k  are  displayed  in  Fig.  5.2.  This  shift  from 

CO 

fk(x)  to  F^(x)  =  /  f^(y)dy  is  immediate  using  an  operational  property  of 
x 

the  Laguerre  transform  [A],  Appendix  C.  The  results  were  compared  with  the 
numerical  results  obtained  from  the  Khinchin-Pollaczek  formula  [7],  again 
via  Laguerre  transform  based  calculation.  The  difference  between  the  two 
results  was  bounded  by  10~  at  k  =  100  and  by  ID”  at  k  =  200. 


Fig.  5.3  displays  similar  results  for  the  G/G/l  queueing  system  with 

_ _ j2 

s(x)  =  /2/n  e  1  ,  the  folded  normal  distribution  (Chi  distribution  with  1 

,  2 

d.f.)  and  t(x)  =  xe  2  ,  the  Rayleigh  distribution.  For  this  G/G/l  setting. 


no  analytical  results  in  the  real  domain  are  available. 


Even  though  the  traffic  intensity  for  the  M/G/l  example  (p  =  0.63)  is 
less  than  that  for  the  G/G/l  example  (p  =  0.64),  the  convergence  in  M/G/l 
is  much  slower.  We  note  that  ?k(0+)  ■+  p  as  k  •+  ®  in  the  M/G/l  example, 
while  this  does  not  hold  in  the  G/G/l  example. 

The  Laguerre  dagger  coefficients  of  the  folded  normal  density  were 
generated  by  an  efficient  recurrence  formula  developed  in  [10] .  Those  for 
the  Rayleigh  density  were  then  derived  using  an  operational  property  of 
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the  Laguerre  transform  [A],  Appendix  C.  In  the  computations  250  terms  were 
used  in  each  of  the  Laguerre  series  (1.1)  for  s(x)  and  t(x). 


00 

/  fk(y)dy 


Fig.  5.2.  The  survival  function  of  the  waiting  time  of  the  k-th  customer 


QO 


Fig.  5.3.  The  survival  function  of  the  waiting  time  of  the  k-th  customer 


Operational  Properties  of  bilateral  Lap.ucrrc  Transforus 


f("0,  g(t)  e  C+(R)  is  assumed.  When  f(-r)  e  C^(R),  the  operational  rules  above 
unchanged  except  for  those  functions  marked  with  *.  The  required  modifications 
given  in  Section  4. 
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APPENDIX  B 

Some  useful  identities  for  the  Laguerre  coefficients  (f*) 

* 

As  before,  f(x)  e  C . (R)  is  assumed. 

▼ 


(1)  f(T)  -  I  fhn(T) 


(2)  fj  =  /  f(T)hn(T)dT 


(3)  If l  *  I  f2(T)dT 


(4]  f#  =  f+  -  f+  ;  f*  «  J  f*  -  -  I 
n  n  n-1  n  L  m 


n+1 


m 


( 5 )  T*(u)  =  l  fV  =  VTtS}  ;  *B(s)  =  /  e~5Tf(T)dx 


.1  1+u. 


-ST, 


(6)  Tj(e16)  =  *g(|  cot  |) 


°°  +  “  # 

(7)  I  f  =  -  l  nfn  -  f (0+)  -  f(O-) 


(8)  I  fn  -  0 


(9)  2  l  (-l)V  »  l  (-l)V  =  /  f(T)dT 


(10)  l  (-l)nnf*  =  j  ]  tf(T)dT 

-00  -CO 

(11)  l  (-DnnV  =  ^/  T2f(T)dx 

-00  -oo 


and  (f*) 
- n— 
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*  (12)  l  f*  =  /  e'^1 T I r f(r)dt  ;  rf(t)  =  f(r)*f(-T) 

_oo  _co  1  1 

*  U3)  )  |n|f*  =  /  | r !f2(r)dT 

-00  —  CO 

*  (14)  l  n2fj  =  i  /  T2f2(T)dT  +  /  T2f'(T)2dT 

—  OO  —  CD  _CO 

*  For  (12),  (13)  and  (14),  f(r)  e  C+(R)  is  assumed. 
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